11institutetext: Harvard Medical School, Boston, MA 02115, USA
11email:  {bian4, awjang<EMAIL_ADDRESS>
# Multi-task Magnetic Resonance Imaging Reconstruction using Meta-learning
††thanks: This work is supported by the National Institute of Health under
grants R21EB031185, R01AR081344, R01AR079442, and R56AR081017.

Wanyu Bian 11 Albert Jang 11 Fang Liu 11

###### Abstract

Using single-task deep learning methods to reconstruct Magnetic Resonance
Imaging (MRI) data acquired with different imaging sequences is inherently
challenging. The trained deep learning model typically lacks generalizability,
and the dissimilarity among image datasets with different types of contrast
leads to suboptimal learning performance. This paper proposes a meta-learning
approach to efficiently learn image features from multiple MR image datasets.
Our algorithm can perform multi-task learning to simultaneously reconstruct MR
images acquired using different imaging sequences with different image
contrasts. The experiment results demonstrate the ability of our new meta-
learning reconstruction method to successfully reconstruct highly-undersampled
k-space data from multiple MRI datasets simultaneously, outperforming other
compelling reconstruction methods previously developed for single-task
learning.

###### Keywords:

Meta-learning MRI Image Reconstruction.

##  1 Introduction

MR images acquired using different imaging sequences can provide complementary
soft tissue contrasts, which are collectively used for clinical disease
diagnosis. However, MRI is a slow imaging modality, resulting in a high
sensitivity to subject motion. Long acquisition time is also undesirable, as
it leads to lower patient throughput than other popular imaging modalities.
Accelerated MRI acquisition and reconstruction are highly desirable and remain
an active research topic [16, 6, 14, 10, 11, 12]. In the past years, deep
learning methods [8] have shown great potential to enable rapid MRI. Deep
learning approaches typically use task-specific deep networks to learn image
features associated with incomplete k-space sampling and aim to remove image
artifacts and noises caused by k-space undersampling during accelerated
acquisition. While the conventional single-task models can work well when the
training and testing data stem from the same data distribution (e.g., acquired
using the same imaging sequence), they typically underperform when training
and testing data differ substantially (e.g., acquired from different imaging
sequences) [1, 15, 9]. Therefore, multi-task deep learning using methods to
synergistically learn image features across multiple image datasets with
different image contrasts is highly desirable for robust and high-efficient
reconstruction.

This paper proposes a meta-learning framework to learn image features from
multiple MR image datasets. Meta-learning is a stacking ensemble learning [4]
scheme that is considered a process of “learning-to-learn”, which can learn to
improve the learning algorithm and parameter generalizability over multiple
training episodes, thus enabling each task to learn better [7]. In this study,
we developed a bi-level meta-learning reconstruction framework (i.e.,
including base-level and meta-level) to handle the learning of multiple image
datasets. At the base-level, we introduce new deep networks (e.g., base-
learners) by unrolling proximal gradient descent in both image and k-space
domains to cross-learn the image and frequency domain features for single
image contrast. In the meta-level, we introduce an optimization algorithm that
can alternatively optimize the base-learners and one additional meta-learner
for efficiently characterizing mutual correlation among multiple image
datasets. Our bi-level meta-learning can simultaneously reconstruct highly-
undersampled k-space data acquired using different imaging sequences with an
optimal reconstruction for all image contrasts.

Our contribution can be summarized as follows:
1\. A novel bi-level meta-learning framework is proposed to reconstruct
highly-undersampled MRI datasets acquired using different imaging sequences.
2\. An unrolled network is proposed to extend proximal gradient descent [13]
on both image and k-space domains, enabling cross-domain learning of MRI data
and leading to a superior reconstruction performance for every single
contrast.
3\. The proposed algorithm is validated for reconstructing a set of knee MRIs
acquired at different imaging contrasts and planes.

##  2 Methodology and Algorithm

Considering the multi-coil MR images at multi-task image datasets, we denote
each image as $\mathbf{x}_{i,j}\in\mathbb{C}^{n}$ for $i=1,\cdots,m$ and
$j=1,\cdots,c$, where $m$ is the number of tasks (i.e., the acquired image
sequences), $c$ is the number of receiving coils and $n$ is the number of
pixels in the image. The corresponding undersampled k-space data can be
formulated as
$\mathbf{f}_{i,j}=\mathbf{P}\mathbf{F}\mathbf{x}_{i,j}+\mathbf{n}_{i,j}$,
where $\mathbf{F}$ denotes the MR encoding matrix, $\mathbf{P}$ is a binary
undersampling matrix, and $\mathbf{n}_{i,j}$ is the noise. We aim to handle
this reconstruction by solving the following optimization problem using meta-
learning:

$\displaystyle\min_{\mathbf{x}_{i,j}}\,\,$
$\displaystyle\Psi_{\Theta,W}(\mathbf{x}_{i,j}):=\textstyle{\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{c}\frac{1}{2}\|\mathbf{P}\mathbf{F}\mathbf{x}_{i,j}-\mathbf{f}_{i,j}\|_{2}^{2}+G_{w_{i}}(\mathbf{x}_{i,j})+K_{w_{i}}(\mathbf{F}(\mathbf{x}_{i,j}))},$
(1a) s.t.
$\displaystyle\mathbf{x}_{i,j}=\operatorname*{arg\,min}_{\mathbf{x}_{i,j}}\textstyle{\frac{1}{2}\|Z_{w_{i}}(H_{\Theta}([J_{w_{i}}(\\{\mathbf{x}_{j}\\})_{i}]_{i=1}^{m}))-RSS(\\{\mathbf{x}_{j}\\})_{i}\|^{2}_{2}}.$
(1b)

In model (1), $G,K,J,Z$ and $H$ all represent deep neural networks. We denote
a collection of task-specific features as $W=\\{w_{i}:i\in[m]\\}$, which are
generated by _base-learner_ s $G,K,J,Z$ to learn from the individual task. The
_meta-learner_ $H$, however, attempts to learn meta-knowledge $\Theta$ that
captures multi-task features from multiple datasets. Note that $W$ and
$\Theta$ are network weights to be learned during training.

More specifically, in the upper-level minimization problem (1a), two learnable
regularization terms, namely $G_{w_{i}}$ and $K_{w_{i}}$, are applied to each
task to learn image domain features and frequency domain (k-space) features,
respectively, from the training data. The learner $G_{w_{i}}$ is designed to
reduce image domain artifacts and noise while the learner $K_{w_{i}}$ intends
to refine k-space to emphasize structural details, patterns, and textures. The
lower-level optimization problem (1b) enforces data consistency to ensure that
the output generated by the three learners $J,H,Z$ are close to the root sum-
of-squares ($RSS$) of the multi-coil images. Each learner performed
fundamentally different but essential roles in the learning process.

$J$ is the coil combination network that integrates the multi-coil images to
produce a coil-combined image [3]. Each task is associated with a task-
specific $J$, and the resulting coil combined image
$J_{w_{i}}(\\{\mathbf{x}_{j}\\})$ for $i$-th task is input into the high-
dimensional meta-learner $H$ to extract meta-knowledge $\Theta$. The meta-
knowledge is generated by learning the mapping
$H_{\Theta}:\mathbb{C}^{mn}\rightarrow\mathbb{C}^{dn}$, where $d$ is the
feature dimension (e.g., the number of kernels in the last convolutional layer
of $H$). The meta-learner $H$ plays an essential role in generating,
balancing, and compensating cross-correlations between the features of the
different tasks. More importantly, the meta-knowledge $\Theta$ provides the
parameter initialization of the entire model, which subsequently supports the
individual base-learners to train their target tasks. The _meta-distributor_
$Z_{w_{i}}:\mathbb{C}^{dn}\rightarrow\mathbb{C}^{n}$ learns the mapping from
the meta-knowledge $\Theta$ to the coil-combined image through $RSS$. More
specifically,
$Z_{w_{i}}(H_{\Theta}([J_{w_{1}}(\\{\mathbf{x}_{1,j}\\}),\cdots,J_{w_{m}}(\\{\mathbf{x}_{m,j}\\})]))$
attempts to output an image matching the coil-combined image
$RSS(\\{\mathbf{x}_{j}\\}_{j=1}^{c})_{i}$ for the $i$-th task. The meta-
distributor maintains data consistency by learning the ability to distribute
the meta-knowledge into each image. This enables knowledge sharing across
multi-tasks, where the meta-knowledge can be efficiently distributed to each
task according to their respective feature importance, thus improving the
reconstruction quality of each task.

In our algorithm, we developed a neural network architecture Fig. 1 to
implement training of the multi-task MRI reconstruction model. Our
implementation consists of a forward step to unroll the optimization procedure
through proximal gradient descent, and a backward meta-training to update
weights and features in base-learners and meta-learner.

Figure 1: The proposed multi-task meta-learning MRI reconstruction.

###  2.1 Forward Learnable Optimization using Gradient Descent

We propose to develop an unrolled network to solve both the upper-level
variational minimization problem and the lower-level constrained optimization
problem in model (1). We introduce a proximal gradient descent [13] inspired
algorithm 1, which can be parametrized as an iterative procedure with total
$T$ iterations. The input $\mathbf{x}^{(0)}_{i,j}$ of the algorithm is
obtained from an initial network that tries to pre-estimate the missing
k-space data (similar to conventional parallel imaging methods such as GRAPPA
[6]) therefore providing a good initial start for Algorithm 1. This helps
reduce the number of iterations needed for convergence and decreases the
overall computational complexity. The overall procedure consists of two for-
loops: an outer and an inner loop, as follows:

1: Input: $\mathbf{x}^{(0)}_{i,j},\delta_{0},\rho_{i}^{(0)}$.

2: for $t=0,1,\dots,T-1$ do

3: for $\tau=0,1,\dots,r$ do

4:
$\mathbf{x}^{(t)}_{i,j}\leftarrow{\mathbf{x}}^{(t)}_{i,j}-\delta_{\tau}\nabla\|Z_{i}(H([J_{i}(\\{{\mathbf{x}}_{j}^{(t)}\\})]_{i\in[m]}))-RSS(\\{\mathbf{x}^{(t)}_{j}\\})_{i}\|$

5: end for

6:
$\mathbf{b}_{i,j}^{(t)}=\mathbf{x}_{i,j}^{(t)}-\rho^{(t)}_{i}\mathbf{F}^{H}\mathbf{P}^{\top}(\mathbf{P}\mathbf{F}\mathbf{x}_{i,j}^{(t)}-\mathbf{f}_{i,j}),\quad
i=1,\cdots,m,\,j=1,\cdots,c,$

7:
$\bar{\mathbf{x}}_{i,j}^{(t)}=[\operatorname{\mathrm{prox}}_{\rho^{(t)}_{i}G_{i}(\cdot)}({\mathbf{b}}^{(t)})]_{i,j},\quad
i=1,\cdots,m,\,j=1,\cdots,c,$

8:
${\mathbf{x}}_{i,j}^{(t+1)}=[\operatorname{\mathrm{prox}}_{\rho^{(t)}_{i}K_{i}(\mathbf{F}(\cdot))}(\bar{\mathbf{x}}^{(t)})]_{i,j},\quad
i=1,\cdots,m,\,j=1,\cdots,c,$

9: end for and output $\mathbf{x}^{(T)}_{i,j}$,
$Z_{i}(H([J_{i}(\\{{\mathbf{x}}_{j}^{(T)}\\})]_{i\in[m]}))$,
$Z_{i}(H([J_{i}(\\{{\mathbf{x}}_{j}^{(0)}\\})]_{i\in[m]}))$

Algorithm 1 Forward Learnable Descent Algorithm for solving (1)

The inner loop aims to solve the lower-level optimization in (1b). We first
calculate the gradient of the constrained least squares in (1b) and then
iterate $r$ steps using gradient descent with a learnable step size
$\delta_{\tau}$ shown in steps 3-5 in Algorithm 1. In the outer loop, the
updated ${\mathbf{x}}^{(t)}_{i,j}$ serves as the initial input for steps 6-8
of proximal gradient descent. The proximal operator is defined as
$\operatorname{\mathrm{prox}}_{\alpha
R}(\mathbf{b}):=\operatorname*{arg\,min}_{\mathbf{x}}\left\\{\frac{1}{2\alpha}\|\mathbf{x}-\mathbf{b}\|^{2}+R(\mathbf{x})\right\\}$.
In step 7 and 8, two proximal operators are conducted by alternating between
image and k-space domains to ensure learning of both image and frequency
features. Because the proximal operators do not have a closed-form solution
for explicit implementation, we parametrize
$\operatorname{\mathrm{prox}}_{\rho^{(t)}_{i}G_{i}(\cdot)}$ and
$\operatorname{\mathrm{prox}}_{\rho^{(t)}_{i}K_{i}(\mathbf{F}(\cdot))}$ as two
learnable deep residual networks in image and k-space domains, respectively,
to emphasize artifacts and noise removal (in image domain) and structure
preservation (in frequency domain). Our final outputs in step 9 from Algorithm
1 are reconstructed multi-coil images of all tasks $\mathbf{x}^{(T)}_{i,j}$,
we then apply $Z_{i},H,J_{i}$ to obtain their corresponding coil-combined,
meta-knowledge distributed image
$\hat{\mathbf{x}}_{i}=Z_{i}(H([J_{i}(\\{{\mathbf{x}}_{j}^{(T)}\\})]_{i\in[m]}))$.
In addition, we also apply $Z_{i},H,J_{i}$ to the initial reconstruction
$\mathbf{x}^{(0)}_{i,j}$ and get
$\hat{\mathbf{x}}_{i}^{(0)}=Z_{i}(H([J_{i}(\\{{\mathbf{x}}_{j}^{(0)})\\}]_{i\in[m]}))$.
Both $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{x}}_{i}^{(0)}$ are then input to
the loss function. Our overall loss function is designed as follows for the
forward optimization:

$\displaystyle\textstyle{\mathcal{L}(\Theta,W;\mathcal{D}):=\sum^{m}_{i=1}\ell_{i}(\Theta,w_{i};\mathcal{D}_{i})},\text{
where each $\ell_{i}$ is called \emph{base-loss}:}$ (2a)
$\displaystyle\textstyle{\ell_{i}(\Theta,w_{i};\mathcal{D}_{i})=\|\hat{\mathbf{x}}_{i}(\Theta,w_{i})-\mathbf{x}^{*}_{i}\|+\lambda\|\hat{\mathbf{x}}_{i}^{(0)}(\Theta,w_{i})-\mathbf{x}^{*}_{i}\|-\mu\text{SSIM}(\hat{\mathbf{x}}_{i}(\Theta,w_{i}),\mathbf{x}^{*}_{i})},$
(2b)

where we use $*$ to reflect the fully-sampled k-space and the ground truth for
each task is $\mathbf{x}^{*}_{i}=RSS(\\{\mathbf{x}^{*}_{j}\\}_{j=1}^{c})_{i}$
and $\mathcal{D}_{i}={(\mathbf{f}_{i},\mathbf{x}^{*}_{i})}$ is the training
data pair for the $i$-th task. The objective is to learn network weights $W$
and $\Theta$ during training. The loss function consists of three terms in
(2). The first term is a standard supervised learning term forcing the final
reconstructed image to approximate the target image from a fully-sampled
k-space. The second term aims to learn a favorable initial input
$\mathbf{x}^{(0)}_{i,j}$ for Algorithm 1. The third term uses the structural
similarity index (SSIM) [17] to enforce the similarity between the
reconstructed and target images. Parameters $\lambda$ and $\mu$ are prescribed
to balance the three terms.

###  2.2 Backward Network Update through Meta-training

The backward operation (3) is conducted during training to minimize the
overall loss function in (2) with respect to the network weights $\Theta$ and
$W$.

$\displaystyle\textstyle{\min_{\Theta(W)}\mathcal{L}(\Theta,W;\mathcal{D})}$
(3a) s.t. $\displaystyle
w_{i}=\textstyle{\operatorname*{arg\,min}_{w_{i}}\ell_{i}(\Theta(W),w_{i};\mathcal{D}_{i}),\quad
i=1,\cdots,m}$ (3b)

The overall procedure to solve (3) is inspired by Model-Agnostic Meta-Learning
(MAML) [5] that consists of two for-loops: an outer and an inner loop, as
follows:

11: Input: $\mathcal{D}_{i}$ and initialize $\Theta$, ${w_{i}}$,
$i=1,\cdots,m$.

2: for $l=0,1,\dots,L$ do

3: Sample training volume $\mathcal{V}_{i}\subset\mathcal{D}_{i}$.

4: for $k=0,1,\dots,K$ (inner loop) do

5:
$\Theta\leftarrow\Theta-\beta\nabla_{\Theta}\mathcal{L}(\Theta,W;\\{\mathcal{V}_{i}:i\in[m]\\})$

6: end for

7: $w_{i}\leftarrow
w_{i}-\alpha_{i}\nabla_{w_{i}}\ell_{i}(\Theta,w_{i};\mathcal{V}_{i}),\quad
i=1,\cdots,m$

8: end for and output $\Theta$, $W=\\{w_{1},\cdots,w_{m}\\}$

Algorithm 2 Backward Meta-Training algorithm for solving (3).

We train the outer loop with a total of $L$ epochs. A mini-batch training with
a randomly selected 3D image volume is used in each epoch
$l\in\\{0,\cdots,L\\}$. The inner loop aims to learn meta-knowledge
$\Theta(W)$ through the upper-level (meta-level) minimization in (3a), where
the meta-knowledge $\Theta$ is a function dependent on network weights $W$
(a.k.a., learning-to-learn). The upper-level minimization fixes $W$ and only
updates $\Theta$ in the inner loop steps 4-6, and then the updated $\Theta(W)$
is used as input into lower-level (base-level) minimization to optimize each
$w_{i}$ in (3b) in the outer loop. $\beta$ represents the meta-learning rate
in the inner loop and $\alpha_{i}$ represents each base learning rate in the
outer loop.

##  3 Experiments

###  3.1 Experiment Setup

The study was approved by the local IRB. Knee datasets with fully-sampled
k-space were acquired on $25$ subjects (20/5 for training/testing). The
datasets were acquired on a 3.0T scanner with an 18-element knee coil array
using four two-dimensional fast spin-echo (FSE) sequences, including coronal
proton density-weighted FSE (Cor-PD) and T2-weighted FSE (Cor-T2), and
sagittal proton density-weighted FSE (Sag-PD) and T2-weighted FSE (Sag-T2).

Experiments were performed to evaluate the efficacy and efficiency of our
proposed method to reconstruct highly-undersampled multi-task knee MRI.
K-space data were retrospectively undersampled [10] to simulate acceleration
rates (AR) of $\\{4,5,6\\}$. To illustrate the effect of meta-learning, we
denote our proposed method without the meta-learner $H_{\Theta}$ (i.e., only
iterate steps 6 to 8 in Algorithm 1 for m tasks separately) as single-task
learning (STL) and the full version as multi-task meta-learning (MTML). All
hyper-parameter selection is shown in the Appendix. We also compared our
proposed STL and MTML against two recently developed reconstruction methods
ISTA-Net [18] and pMRI-Net [2], both of which were previously shown to provide
great performance in reconstructing undersampled single-task MRI data.

###  3.2 Experimental Results

Table 1: Quantitative comparison of different methods. (ITL: proposed single-
task learning; MTML: proposed multi-task meta-learning; PSNR: peak signal-to-
noise ratio; NMSE: normalized mean squared error)

|  | Sag-T2 |
---|---|---|---
AR | 4x | 5x | 6x
| PSNR/SSIM/NMSE | PSNR/SSIM/NMSE | PSNR/SSIM/NMSE
ISTA-Net[18] | 34.6138/0.9687/0.1444 | 33.4308/0.9600/0.1659 | 31.5700/0.9458/0.2061
pMRI-Net[2] | 36.1084/0.9769/0.1214 | 35.0315/0.9712/0.1375 | 34.3608/0.9672/0.1486
ITL | 37.4666/0.9825/0.1039 | 36.6242/0.9793/0.1144 | 35.9206/0.9763/0.1239
MTML | 38.3446/0.9854/0.0941 | 37.5966/0.9832/0.1023 | 36.8210/0.9799/0.1118

|  | Cor-T2 |
---|---|---|---
| PSNR/SSIM/NMSE | PSNR/SSIM/NMSE | PSNR/SSIM/NMSE
ISTA-Net[18] | 34.1552/0.9628/0.1142 | 32.7243/0.9472/0.1382 | 30.0678/0.9155/0.1864
pMRI-Net[2] | 36.3567/0.9767/0.0867 | 35.5938/0.9728/0.0942 | 34.1399/0.9646/0.1108
ITL | 37.0114/0.9796/0.0808 | 36.1727/0.9758/0.0883 | 35.7136/0.9739/0.0930
MTML | 37.8478/0.9828/0.0740 | 37.1512/0.9812/0.0793 | 36.6107/0.9781/0.0843

|  | Sag-PD |
---|---|---|---
| PSNR/SSIM/NMSE | PSNR/SSIM/NMSE | PSNR/SSIM/NMSE
ISTA-Net[18] | 33.1395/0.9541/0.0774 | 31.9699/0.9495/0.0877 | 29.4236/0.9280/0.1176
pMRI-Net[2] | 35.0189/0.9685/0.0623 | 34.1714/0.9588/0.0690 | 31.8230/0.9365/0.0906
ITL | 35.9640/0.9736/0.0558 | 35.0720/0.9680/0.0617 | 34.0919/0.9619/0.0691
MTML | 37.0246/0.9784/0.0493 | 36.2652/0.9754/0.0537 | 35.1452/0.9688/0.0611

|  | Cor-PD |
---|---|---|---
| PSNR/SSIM/NMSE | PSNR/SSIM/NMSE | PSNR/SSIM/NMSE
ISTA-Net[18] | 32.3933/0.9640/0.1030 | 31.1878/0.9555/0.1150 | 29.1795/0.9603/0.1281
pMRI-Net[2] | 33.5534/0.9695/0.0937 | 32.7408/0.9667/0.1003 | 30.8576/0.9560/0.1178
ITL | 35.5481/0.9760/0.0809 | 33.5579/0.9698/0.0939 | 31.9874/0.9625/0.1045
MTML | 36.6621/0.9797/0.0748 | 35.4338/0.9765/0.0817 | 32.8284/0.9663/0.0999

Figure 2: Comparison of Sag-T2 and Cor-T2 contrasts at AR = 6.

Fig. 2 compares the reconstructed images for Sag-T2 and Cor-T2 contrasts from
different methods at AR=6. While ISTA-Net and pMRI-Net can largely remove the
aliasing artifacts at this undersampling level, both methods yielded
suboptimal reconstruction performance with noticeable image blurring in the
reconstructed images. Our proposed STL can enable better reconstruction than
ISTA-Net and pMRI-Net with improved image sharpness, due to its ability to
cross-learn both image and k-space features, therefore better removing noises
and artifacts and meanwhile preserving high-frequency image features. Our
proposed MTML using meta-learning performs the best, outperforming all single-
task learning methods. The pixel-wise error maps show that MTML has the
smallest reconstruction errors compared to the fully-sampled reference, and
the error is more homogeneous across the entire image domain for MTML. The
zoom-in images further highlight the superb reconstruction performance of MTML
where the tissue boundaries are better characterized at all tissue types
(cartilage, meniscus, and muscle), and tissue texture, sharpness, and
conspicuity are well-preserved. Similar comparisons of reconstructed images
for the Sag-PD and Cor-PD sequences are shown in the Appendix, where MTML also
consistently enables the best and most balanced performance on different types
of contrast.

The qualitative comparison in Fig. 2 is further verified by different
quantitative metrics shown in Table 1. The table shows the mean PSNR, SSIM,
and NMSE for all tested subjects for all methods at three ARs=$\\{4,5,6\\}$.
MTML achieves the best quantitative performance in all metrics at all ARs,
consistent with the qualitative assessment in exemplified figures.

Our proof-of-concept study using multi-task meta-learning opens a new window
for further investigating robust and efficient multi-contrast large-scale MRI
reconstruction algorithms.

## References

  * [1] Antun, V., Renna, F., Poon, C., Adcock, B., Hansen, A.C.: On instabilities of deep learning in image reconstruction and the potential costs of ai. Proceedings of the National Academy of Sciences 117(48), 30088–30095 (2020)
  * [2] Bian, W., Chen, Y., Ye, X.: Deep parallel mri reconstruction network without coil sensitivities. In: Machine Learning for Medical Image Reconstruction: Third International Workshop, MLMIR 2020, Held in Conjunction with MICCAI 2020, Lima, Peru, October 8, 2020, Proceedings 3. pp. 17–26. Springer (2020)
  * [3] Bian, W., Chen, Y., Ye, X.: An optimal control framework for joint-channel parallel mri reconstruction without coil sensitivities. Magnetic Resonance Imaging (2022)
  * [4] Džeroski, S., Ženko, B.: Is combining classifiers with stacking better than selecting the best one? Machine learning 54, 255–273 (2004)
  * [5] Finn, C., Abbeel, P., Levine, S.: Model-agnostic meta-learning for fast adaptation of deep networks. In: International conference on machine learning. pp. 1126–1135. PMLR (2017)
  * [6] Griswold, M.A., et al.: Generalized autocalibrating partially parallel acquisitions (grappa). Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 47(6), 1202–1210 (2002)
  * [7] Hospedales, T., Antoniou, A., Micaelli, P., Storkey, A.: Meta-learning in neural networks: A survey. IEEE transactions on pattern analysis and machine intelligence 44(9), 5149–5169 (2021)
  * [8] Knoll, F., Hammernik, K., Zhang, C., Moeller, S., Pock, T., Sodickson, D.K., Akcakaya, M.: Deep-learning methods for parallel magnetic resonance imaging reconstruction: A survey of the current approaches, trends, and issues. IEEE signal processing magazine 37(1), 128–140 (2020)
  * [9] Liu, F., Samsonov, A., Chen, L., Kijowski, R., Feng, L.: Santis: sampling-augmented neural network with incoherent structure for mr image reconstruction. Magnetic resonance in medicine 82(5), 1890–1904 (2019)
  * [10] Lustig, M., Donoho, D., Pauly, J.M.: Sparse mri: The application of compressed sensing for rapid mr imaging. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 58(6), 1182–1195 (2007)
  * [11] Lustig, M., Pauly, J.M.: Spirit: iterative self-consistent parallel imaging reconstruction from arbitrary k-space. Magnetic resonance in medicine 64(2), 457–471 (2010)
  * [12] Otazo, R., Kim, D., Axel, L., Sodickson, D.K.: Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion mri. Magnetic resonance in medicine 64(3), 767–776 (2010)
  * [13] Parikh, N., Boyd, S.: Proximal algorithms. Foundations and Trends in optimization 1(3), 127–239 (2014)
  * [14] Pruessmann, K.P., Weiger, M., Scheidegger, M.B., Boesiger, P.: Sense: sensitivity encoding for fast mri. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine 42(5), 952–962 (1999)
  * [15] Shimron, E., Tamir, J.I., Wang, K., Lustig, M.: Implicit data crimes: Machine learning bias arising from misuse of public data. Proceedings of the National Academy of Sciences 119(13), e2117203119 (2022)
  * [16] Sodickson, D.K., Manning, W.J.: Simultaneous acquisition of spatial harmonics (smash): fast imaging with radiofrequency coil arrays. Magnetic resonance in medicine 38(4), 591–603 (1997)
  * [17] Wang, Z., et al.: Image quality assessment: from error visibility to structural similarity. IEEE transactions on image processing 13(4), 600–612 (2004)
  * [18] Zhang, J., Ghanem, B.: Ista-net: Interpretable optimization-inspired deep network for image compressive sensing. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 1828–1837 (2018)


#  Frictional contact of soft polymeric shells

Riad Sahli  0000-0001-9077-3602 Department of Mechanical and Production
Engineering, Aarhus University, Denmark  Jeppe Mikkelsen  Department of
Mechanical and Production Engineering, Aarhus University, Denmark  Mathias
Sätherström Boye  Department of Mechanical and Production Engineering, Aarhus
University, Denmark  Marcelo A. Dias  0000-0002-1668-0501
<EMAIL_ADDRESS>Institute for Infrastructure and Environment, School of
Engineering, The University of Edinburgh, EH9 3FG Edinburgh, UK  Ramin
Aghababaei  0000-0002-0700-0084<EMAIL_ADDRESS>Department of Mechanical and
Production Engineering, Aarhus University, Denmark

###### Abstract

The classical Hertzian contact model establishes a monotonic correlation
between contact force and area. Here, we showed that the interplay between
local friction and structural instability can deliberately lead to
unconventional contact behavior when a soft elastic shell comes into contact
with a flat surface. The deviation from Hertzian contact first arises from
bending within the contact area, followed by the second transition induced by
buckling, resulting in a notable decrease in the contact area despite
increased contact force. Friction delays both transitions and introduces
hysteresis during unloading. However, a high amount of friction suppresses
both buckling and dissipation. Different contact regimes are discussed in
terms of rolling and sliding mechanisms, providing insights for tailoring
contact behaviors in soft shells.

a) b) |  c)
---|---
Figure 1: (a) Bulk hemisphere made of PDMS in contact with a PLA plate. The
contact area is circular and grows with normal load with an exponent of 2/3,
following Hertzian contact solution (red line in subfig. c). (b) PDMS
hemispherical shell in dry contact with a PLA plate (solid blue line in
subfig. c.). (c) Presents the variation of contact area ($\bar{A}=2A/(\pi
R^{2})$) as a function of normal load ($\bar{F}=FR/(Eh^{3})$) for bulk and
shell. First, a smooth deviation from Hertzian contact can be observed
(inset). As indentation progresses, the contact morphology transitions from
circular to annular, causing a sudden drop in the contact area (marked by a
double arrow).

The contact area and pressure distribution between solid surfaces directly
influence the frictional [1, 2] and thermal resistance [3] as well as the
sealing capacity [4] of the interface. Assuming a parabolic pressure profile,
infinitesimal deformation, and material homogeneity, Hertzian theory [5]
establishes a power-law relationship between the normal load and contact area
for a wide range of elastic geometries, including spherical, cylindrical, and
elliptical contacts [6] (also see Fig. 1). Alternatively, it is shown that the
contact area increases linearly with the normal load when plastic deformation
occurs [1]. Since these findings, the contact area-load relation faced a
transformative modification as the interplay of other factors such as friction
[7, 8, 9], adhesion [10, 11, 12], roughness [13, 14, 15, 16], and material and
geometric nonlinearities [17, 18, 19] were examined. It is well established
that the contact area grows with the normal load [20], while material
properties, contact geometry, and surface conditions dictate the growth rate
[21, 22, 23, 24].

Recent studies [25, 26] have discussed the possibilities of modifying friction
laws (i.e. the correlation between the normal and frictional forces) by
designing metainterfaces with certain roughness morphology. Studying the
contact response of soft polymeric shells has offered another angle to modify
the contact laws using the structural instability at the contact. Pauchard and
Rica demonstrated a sub-critical bifurcation, namely a first-order phase
transition, changes the nature of the mechanical response in contact [27]. It
was observed that a flat contact exists initially between the rigid plane and
the shell. However, it transitions to a post-buckling state characterized by
contact along an axisymmetric circular ridge formed between the rigid plane
and the inverted shell. The relevant parameters of this problem are material,
shell’s Young’s modulus $E$, and geometrical, shell’s radius $R$ and thickness
$h$. The combination of these results into the relevant force scale,
$F_{0}\equiv Eh^{3}/R$, and a detailed analysis proposes three regimes [28]:
(i) Hertzian regime with a circilar contact area, for a range of forces $F\ll
F_{0}$, where deformation occurs within a region smaller than the shell
thickness, $h$; (ii) intermediate regime, where $F\sim F_{0}$ and equilibrium
is obtained by a balance between bending and stretching, the shell flattens in
a manner assumed to be similar to Hertzian contact with a disk-like contact
area; (iii) post-buckling regimes, where $F>F_{0}$ and the stretching energy
caused by sphere flattening becomes significant, destabilizing the shell and
buckling into an inverted shape, which is akin to the isometric problem
presented by Pogorelov [29]. Recent studies [30, 31, 32] have demonstrated
secondary buckling instability, where a transition from axisymmetric to
asymmetric deformation mode with multiple vertices occurs. Additionally, the
effect of structural imperfections on the buckling mode has also been
investigated [33]. Despite these advances in understanding the mechanics of
shell indentation, there still exist limitations in the exploration of the
effect of structural instability and friction on contact morphology and
pressure distribution.

In this letter, we examine the evolution of contact morphology and pressure
when a thin shell comes into contact with a rigid, flat surface at different
stages of indentation with axisymmetric deformation. Our experiments and
numerical simulations show that a circular contact area is first established
and grows by indentation, followed by a smooth transition from circular to
disk-like contact as a result of bending in the contact area. As indentation
progresses, a secondary transition occurs, resulting in a sharp reduction in
the contact area despite increased contact force. Our results led to three
novel discoveries: (1) the force-contact area deviates from Hertzian solution
in the intermediate regime, where the contact geometry transitions from
circular to disk-like contact; (2) a non-monotonic force-contact area relation
can be achieved in connection with the shell instability (post-buckling
regime), which is dictated by the local friction; and (3) friction introduces
hysteresis upon unloading with a maximum at a critical coefficient of friction
(CoF), after which hysteresis drops as high friction suppresses sliding. One
critical aspect of our findings, which pertains to the contact behavior in the
intermediate regime is that we observed a departure from previous claims that
full contact is maintained until the bifurcation point [28]. Additionally, it
is also shown that local friction delays and eventually suppresses the post-
bifurcation regime, hindering the abrupt change in the force-contact area
relation and hysteresis.

Figure 2: The evolution of (a) normalized force $\bar{F}$, (b) contact area
$\bar{A}$, and (c) maximum contact pressure multiplied by $10^{3}$ for
clarity, ($\bar{p}=P_{max}/E$) as a function of indentation ($\bar{d}=d/R$)
for different CoF. Black diamond and red triangle symbols in (a) present,
respectively, the experimental data under dry and lubricated conditions. Solid
lines represent simulation results for $\mu=0$ (blue), $0.48$ (red), $1.0$
(black), and $2.2$ (green). The inset in (a) presents the dissipated energy
during unloading $W$, highlighting a parabolic correlation with CoF. The inset
in (b) reveals two transitional points, 1) the deviation of force-area
relation from the Hertzian solution (the orange dash line) due to bending, and
2) a sudden drop in contact area due to shell buckling.

The snapshots in Fig. 1a and b showcase the deformation modes of a bulk solid
and a shell, both made from cross-linked polydimethylsiloxane (PDMS), upon
contact with a flat polylactic acid (PLA) plate. Fig. 1c presents the force-
contact area relation, highlighting an early deviation from Hertzian contact,
followed by a sudden drop in the contact area when a transition occurs in the
deformation mode. First, we experimentally study the contact of a thin,
elastic hemispherical shell made of cross-linked PDMS [34] with a flat PLA
plate, having a roughness of $Sa=15$ $\mu m$ to minimize the adhesion effect
(Fig. 1). The PDMS shell, with a radius of $R=25mm$ and thickness $h=1.33\pm
0.06$ mm, is made by pouring PDMS mixture over a $25$ mm radius steel ball
[35]. In a displacement-control manner (see supplementary Fig.
LABEL:figSI:fullsetup [36]), the PLA plate comes into contact with the shell
at a constant velocity of $10$ mm/min, where the force-displacement data are
recorded. A side camera is used to monitor the deformation of the shell. To
study the effect of friction coefficient, we ran the test under dry and
lubricated conditions, where “Kema GL-68” is used as the lubricant. Additional
experiments are conducted to measure $\mu$ between PDMS and PLA under dry
($1.0\pm 0.2$) and lubricated ($0.48\pm 0.09$) contacts (see supplementary
Fig. LABEL:figSI:friction [36]).

a)

b)

c)

Figure 3: (a) In-situ snapshots of shell profile indented by a PLA plate,
under dry contact ($\mu=1.0$) at different stages of indentation,
corresponding to points $p_{0}$-$p_{5}$, marked in Fig. 2a. The profiles
obtained from the simulation are superimposed with red solid lines for
comparison. (b) and (c) present, respectively, the corresponding contact
morphology and the contact pressure distribution ($p(x)/E$), obtained from
simulations multiplied by $10^{3}$ for clarity. Inset numbers in (b) present
the total contact area. The evolution of contact area and pressure
distribution, as well as the indentation force and shell profile for different
CoFs can be seen in supplementary movies [36]).

To complement our experiments with details of contact area and pressure
evolution over a wider range of $\mu$, we performed systematic finite element
simulations using the commercial software “Abaqus”[37]. The numerical setup,
shown in supplementary Fig. LABEL:figSI:abacus [36], simulated the contact
between an axisymmetric PDMS shell with a rigid plate. Details of loading and
boundary conditions as well as material properties are presented in the
Supplementary section [36].

Fig. 2a summarizes the experimental and numerical results, presenting the
evolution of contact force ($\bar{F}=FR/(Eh^{3})$) during the loading and
unloading stages under dry ($\mu=1.0$) and lubricated ($\mu=0.48$) conditions.
We additionally include simulation results for the frictionless ($\mu=0.0$)
and high-friction ($\mu=2.2$) scenarios. Despite the absence of the adhesion
effect in our simulation, a close agreement between the experimental and
numerical results has been obtained for the dry and lubricated contacts in all
stages of indentation. As seen, the contact force monotonically increases with
the indentation until buckling occurs, after which a small drop is identified
in the indentation force. The force continues to increase at a lower rate in
the post-buckling regime. One can see that the evolution of force in the post-
buckling regime is a function of CoF. As friction rises, it opposes the
compressive stresses within the shell, increasing the critical force needed
for instability initiation. In other words, the heightened frictional forces
hinder relative motion between surfaces, strengthening resistance to
compressive loading. Consequently, the critical force threshold rises due to
the combined effects of friction-induced resistance and the reaction to
compressive stresses on the shell material. This observation highlights how
frictional mechanics control the structural stability of materials.

A significant dissipation through hysteresis is also observed upon unloading
in both cases. It is evident that as CoF increases, the buckling is delayed,
while a larger energy dissipation occurs upon unloading. As expected, no
dissipation occurs in the frictionless case, indicating that sliding in the
post-buckling stage is the main source of energy dissipation. Interestingly,
it can be seen that high CoF can completely avoid shell buckling and,
consequently, energy dissipation is reduced to a minimum. The inset of Fig. 2a
show a parabolic relation between the dissipated energy ($W=\oint_{S}Fdx$) and
CoF, indicating the existence of a critical CoF, where sliding and dissipation
in the post-buckling stage are maximized (also see Table.LABEL:tab:1 [36]).
This critical CoF may be a function of geometrical and material parameters.

The evolution of contact area in different stages of indentation is presented
in Fig. 2b. It can be seen that the normalized contact area ($\bar{A}=2A/(\pi
R^{2})$) increases monotonically until the point of instability, after which a
sharp reduction occurs (also see Table.LABEL:tab:1 [36]. For the case of high
friction ($\mu=2.2$) and the absence of instability, the contact area
increases for the entire simulation. The inset figure in Fig. 2b, compares the
force-contact area relationship with the Hertzian solution shown by an orange
dashline in the loading regime. Interestingly, it can be seen that deviation
from the Hertzian solution occurs in the intermediate regime and before the
point of buckling. This observation challenges the previous assumption of
Hertzian-like behavior in the intermediate regime [28]. Remarkably, the
contact area remains constant in the post-buckling regime despite the increase
in the indentation force. This stands in contrast to the contact response of a
bulk sphere, where the contact area monotonically grows [5, 21] (also see Fig.
1). Three regimes of deformation can also be distinguished by monitoring the
average contact pressure as a function of indentation. Fig. 2c presents the
maximum contact pressure normalized by elastic modulus ($\bar{p}=P_{avg}/E$).
The pressure initially increases in the Hertzian regime (regime I), followed
by a smooth drop and saturation (regime II). Clearly, it can be seen that the
pressure in the intermediate regime II does not follow the Hertzian behavior
shown in regime I. The maximum pressure experiences a sharp increase as a
result of buckling and the reduction in the contact area (see Fig. 2b),
indicating the transition into the post-buckling regime III.

An in-situ analysis of shell deformation and the evolution of contact area and
pressure for the dry contact scenario ($\mu=1.0$), corresponding to marked
points $p_{0}$-$p_{5}$ in Fig. 2a, is illustrated in Fig. 3 (see supplementary
movies [36] for other CoFs). It can be seen that the numerical simulations
predict the profile of the shell in all regimes. At the onset of contact
($p_{0}$), the circular contact geometry (Fig. 3b) and the parabolic contact
pressure distribution (Fig. 3c) follow the Hertzian solution. $p_{1}$
demonstrates the deviation from the Hertzian solution in the intermediate
regime but before the point of buckling by showing a disk-like contact area
and non-parabolic pressure distribution as a direct result of bending in the
shell. A similar pressure profile has also been derived theoretically [38].
Contrary to the parabolic Hertzian pressure distribution, the pressure in the
intermediate regime reduces from the leading edge to the trailing edge (see
Supp. Movies [36] ), indicating that the change in the contact morphology is
dominated by the rolling mechanism. This explains the minimization of
dissipation until the point of buckling, after which sliding and thus
dissipation occur at the contact. In the post-buckling regime ($p_{2}$ and
$p_{3}$), the contact occurs along a circular ridge with a parabolic pressure
distribution across the contact width. In the unloading phase, the outer side
of the shell is first unloaded elastically ($p_{4}$) as shown in Supp. Fig.
LABEL:figSI:Fig2. In this regime, while the buckled (inner) side of the shell
remains fairly stationary due to friction (see Supp. Fig.
LABEL:figSI:radiusa), the outer side unrolls, causing a reduction in the
contact area (Fig. LABEL:figSI:radiusb). The unloading phase continues with a
transition, where the local tangential load overcomes friction. In this regime
($p_{5}$), while the contact area remains constant (see Fig. 2b), dissipation
occurs via frictional sliding. Eventually, the shell unbuckles, after which
the Hertzian solution is recovered.

Figure 4: The evolution of contact size and location in the loading and
unloading phases for different CoF. The indentation on the y axes is
normalized by the maximum indentation. The contact size represents the radius
of a circular contact in the Hertzian regime and the width of an annular
contact area in the intermediate and post-buckling regimes. Letters H, I, and
B denote Hertzian, intermediate, and post-buckling regimes respectively. Red
dash lines mark the transition point between the Hertzian and intermediate
regimes. The shaded areas highlight the elastic unloading at the onset of
unloading for cases of $\mu=0.48$ and $1.0$.

Another intriguing effect of friction is that it makes the loading-unloading
response asymmetric. Fig. 4 demonstrates the evolution of contact morphology
during the loading and unloading phases by showing the location and size of
the contact area. The contact size refers to the radius of a circular contact
in the Hertzian regime and the width of the annular contact area in the
intermediate and post-buckling regimes. Noting the variation of contact size,
three regimes of Hertzian, intermediate, and post-buckling are classified for
each friction case. It can be seen that friction postpones all transition
points during the loading and unloading phases, resulting in a bigger Hertzian
regime during the loading, but a smaller one upon unloading. Additionally, the
degree of asymmetry in the loading-unloading response increases with friction.
For the high friction case ($\mu=2.2$) where the buckling regime is absent,
unloading is mainly accommodated via the rolling mechanism, explaining the
absence of dissipation (also see Supp. Movie [36]).

We have shown how the interplay between structural instability and the local
friction coefficient dictates the morphology of contact. It is shown that the
contact behavior departs from the Hertzian solution due to bending at the
onset of contact. This observation challenges the previous assumption of a
Hertzian-like response in the pre-buckling regime. In agreement with previous
studies [28], two transition points are identified in the loading phase,
classifying three regimes of deformation: Hertzian, intermediate, and post-
buckling. While the contact area increases monotonically in the first two
regimes, it abruptly drops and then remains constant in the post-buckling
regime. The distribution of pressure in the contact zone revealed two distinct
mechanisms for the evolution of the contact area: (i) the energy-preserving
rolling mechanism in the Hertzian and intermediate regimes, and (ii) the
dissipative frictional sliding in the post-buckling regime. In other words,
the response is fully reversible (non-dissipative) in the absence of friction
or the post-buckling regime. The unloading phase for $\mu=0.48$ begins with an
elastic unrolling, followed by dissipative frictional sliding. Furthermore, it
is evident that friction causes a delay at all transition points, leading to
the expansion and contraction of the Hertzian regime during the loading and
unloading phases, respectively.

The interplay between friction, material properties, and structural
instability leverages the controlled departure from Hertzian behavior and
monotonic contact force-area relation. This new understanding presents a
direction to tailor the contact morphology and governing mechanism, offering
solutions to a diverse range of engineering applications such as systems with
superior damping capabilities [39], sport ball [40, 41], switch haptic [42,
43], or vesicle manipulation [44, 45].

Acknowledgments. This research has been supported by a Marie Sklodowska-Curie
Postdoctoral Fellowship under contract number (proposal number 101065669).
R.A. and R.S. also acknowledge the support from the Thomas B. Thriges Fond.
M.A.D. would like to thank UKRI for support under the EPSRC Open Fellowship
scheme (Project No. EP/W019450/1).

## References

  * Bowden  _et al._ [1950] F. P. Bowden, D. Tabor, and F. Palmer,  _The friction and lubrication of solids_ (Clarendon Press, Oxford, 1950).
  * Archard [1957] J. Archard, Elastic deformation and the laws of friction, Proceedings of the royal society of London. Series A. Mathematical and physical sciences 243, 190 (1957).
  * Madhusudana and Madhusudana [1996] C. V. Madhusudana and C. Madhusudana,  _Thermal contact conductance_ , Vol. 79 (Springer, 1996).
  * Persson  _et al._ [2004] B. N. Persson, O. Albohr, U. Tartaglino, A. Volokitin, and E. Tosatti, On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion, Journal of physics: Condensed matter 17, R1 (2004).
  * Hertz [1882] H. Hertz, Ueber die berührung fester elastischer körper., (1882).
  * Popov  _et al._ [2019] V. L. Popov, M. Heß, and E. Willert,  _Handbook of contact mechanics: exact solutions of axisymmetric contact problems_ (Springer Nature, 2019).
  * Carlsson  _et al._ [2000] S. Carlsson, S. Biwa, and P.-L. Larsson, On frictional effects at inelastic contact between spherical bodies, International Journal of Mechanical Sciences 42, 107 (2000).
  * Storåkers and Elaguine [2005] B. Storåkers and D. Elaguine, Hertz contact at finite friction and arbitrary profiles, Journal of the Mechanics and Physics of Solids 53, 1422 (2005).
  * Weber  _et al._ [2022] B. Weber, J. Scheibert, M. P. de Boer, and A. Dhinojwala, Experimental insights into adhesion and friction between nominally dry rough surfaces, MRS Bulletin 47, 1237 (2022).
  * Pastewka and Robbins [2014] L. Pastewka and M. O. Robbins, Contact between rough surfaces and a criterion for macroscopic adhesion, Proceedings of the National Academy of Sciences 111, 3298 (2014).
  * Bense  _et al._ [2020] H. Bense, M. Tani, M. Saint-Jean, E. Reyssat, B. Roman, and J. Bico, Elastocapillary adhesion of a soft cap on a rigid sphere, Soft Matter 16, 1961 (2020).
  * Oliver  _et al._ [2023] C. Oliver, D. Dalmas, and J. Scheibert, Adhesion in soft contacts is minimum beyond a critical shear displacement, Journal of the Mechanics and Physics of Solids 181, 105445 (2023).
  * Greenwood and Tripp [1967] J. A. Greenwood and J. H. Tripp, The Elastic Contact of Rough Spheres, Journal of Applied Mechanics 34, 153 (1967).
  * Luan and Robbins [2005] B. Luan and M. O. Robbins, The breakdown of continuum models for mechanical contacts, Nature 435, 929 (2005).
  * Aghababaei  _et al._ [2022] R. Aghababaei, E. E. Brodsky, J.-F. Molinari, and S. Chandrasekar, How roughness emerges on natural and engineered surfaces, MRS Bulletin 47, 1229 (2022).
  * Müser and Nicola [2022] M. H. Müser and L. Nicola, Modeling the surface topography dependence of friction, adhesion, and contact compliance, MRS Bulletin 47 (2022).
  * Weber  _et al._ [2018] B. Weber, T. Suhina, T. Junge, L. Pastewka, A. Brouwer, and D. Bonn, Molecular probes reveal deviations from amontons’ law in multi-asperity frictional contacts, Nature communications 9, 888 (2018).
  * Wu  _et al._ [2016] C.-E. Wu, K.-H. Lin, and J.-Y. Juang, Hertzian load–displacement relation holds for spherical indentation on soft elastic solids undergoing large deformations, Tribology International 97, 71 (2016).
  * Ghaednia  _et al._ [2017] H. Ghaednia, X. Wang, S. Saha, Y. Xu, A. Sharma, and R. L. Jackson, A Review of Elastic–Plastic Contact Mechanics, Applied Mechanics Reviews 69, 060804 (2017).
  * Dieterich and Kilgore [1994] J. H. Dieterich and B. D. Kilgore, Direct observation of frictional contacts: New insights for state-dependent properties, Pure and applied geophysics 143, 283 (1994).
  * Sahli  _et al._ [2018] R. Sahli, G. Pallares, C. Ducottet, I. Ben Ali, S. Al Akhrass, M. Guibert, and J. Scheibert, Evolution of real contact area under shear and the value of static friction of soft materials, Proceedings of the National Academy of Sciences 115, 471 (2018).
  * Sahli  _et al._ [2019] R. Sahli, G. Pallares, A. Papangelo, M. Ciavarella, C. Ducottet, N. Ponthus, and J. Scheibert, Shear-induced anisotropy in rough elastomer contact, Physical Review Letters 122, 214301 (2019).
  * Matsuda  _et al._ [2018] K. Matsuda, S. Mori, A. Hatanaka, T. Sunahara, and K. Nakamura, Effect of specimen thickness on growth of real contact area of rubber with two-dimensional regular wavy surface, Tribology International 124, 184 (2018).
  * Mo  _et al._ [2009] Y. Mo, K. Turner, and I. Szlufarska, Friction laws at the nanoscale, Nature 457, 1116 (2009).
  * Slesarenko and Pastewka [2024] V. Slesarenko and L. Pastewka, The bumpy road to friction control, Science 383, 150 (2024).
  * Aymard  _et al._ [2024] A. Aymard, E. Delplanque, D. Dalmas, and J. Scheibert, Designing metainterfaces with specified friction laws, Science 383, 200 (2024).
  * Pauchard and Rica [1998] L. Pauchard and S. Rica, Contact and compression of elastic spherical shells: the physics of a ‘ping-pong’ball, Philosophical Magazine B 78, 225 (1998).
  * Audoly and Pomeau [2010] B. Audoly and Y. Pomeau,  _Elasticity and geometry: from hair curls to the non-linear response of shells_ (Oxford university press, 2010).
  * Pogorelov [1988] A. V. Pogorelov,  _Bendings of surfaces and stability of shells_ , Vol. 72 (American Mathematical Soc., 1988).
  * Vaziri [2009] A. Vaziri, Mechanics of highly deformed elastic shells, Thin-Walled Structures 47, 692 (2009).
  * Nasto  _et al._ [2013] A. Nasto, A. Ajdari, A. Lazarus, A. Vaziri, and P. M. Reis, Localization of deformation in thin shells under indentation, Soft Matter 9, 6796 (2013).
  * Knoche and Kierfeld [2014] S. Knoche and J. Kierfeld, The secondary buckling transition: wrinkling of buckled spherical shells, The European Physical Journal E 37, 1 (2014).
  * Lee  _et al._ [2016a] A. Lee, F. López Jiménez, J. Marthelot, J. W. Hutchinson, and P. M. Reis, The Geometric Role of Precisely Engineered Imperfections on the Critical Buckling Load of Spherical Elastic Shells, Journal of Applied Mechanics 83, 111005 (2016a).
  * Le Rouzic  _et al._ [2009] J. Le Rouzic, P. Delobelle, P. Vairac, and B. Cretin, Comparison of three different scales techniques for the dynamic mechanical characterization of two polymers (pdms and su8), The European Physical Journal-Applied Physics 48 (2009).
  * Lee  _et al._ [2016b] A. Lee, P.-T. Brun, J. Marthelot, G. Balestra, F. Gallaire, and P. M. Reis, Fabrication of slender elastic shells by the coating of curved surfaces, Nature communications 7, 1 (2016b).
  * [36] Supplementary information.
  * Abaqus [2021] G. Abaqus, Abaqus 15.50, Dassault Systemes Simulia Corporation, Providence, RI, USA , 3 (2021).
  * Updike and Kalnins [1972] D. P. Updike and A. Kalnins, Contact pressure between an elastic spherical shell and a rigid plate, Journal of Applied Mechanics 39 (1972).
  * Ullah  _et al._ [2022] R. Ullah, M. Vafaei, S. C. Alih, and A. Waheed, A review of buckling-restrained braced frames for seismic protection of structures, Physics and Chemistry of the Earth, Parts A/B/C , 103203 (2022).
  * Rémond  _et al._ [2022] T. Rémond, V. Dolique, F. Vittoz, S. Antony, R. G. Rinaldi, L. Manin, and J.-C. Géminard, Dynamical buckling of a table-tennis ball impinging normally on a rigid target: Experimental and numerical studies, Physical Review E 106, 014207 (2022).
  * Rinaldi  _et al._ [2016] R. G. Rinaldi, L. Manin, C. Bonnard, A. Drillon, H. Lourenco, and N. Havard, Non linearity of the ball/rubber impact in table tennis: experiments and modeling, Procedia engineering 147, 348 (2016).
  * Plotnick [2017] R. Plotnick, Force, flatness and touch without feeling: Thinking historically about haptics and buttons, New Media & Society 19, 1632 (2017).
  * Gedsun  _et al._ [2022] A. Gedsun, R. Sahli, X. Meng, R. Hensel, and R. Bennewitz, Bending as key mechanism in the tactile perception of fibrillar surfaces, Advanced Materials Interfaces 9, 2101380 (2022).
  * Quemeneur  _et al._ [2012] F. Quemeneur, C. Quilliet, M. Faivre, A. Viallat, and B. Pépin-Donat, Gel phase vesicles buckle into specific shapes, Physical review letters 108, 108303 (2012).
  * Tsapis  _et al._ [2005] N. Tsapis, E. R. Dufresne, S. S. Sinha, C. S. Riera, J. W. Hutchinson, L. Mahadevan, and D. A. Weitz, Onset of buckling in drying droplets of colloidal suspensions, Physical review letters 94, 018302 (2005).


# Pole properties of a resonance: When to subtract partial-decay widths to
obtain the pole widths††thanks: Presented at Excited QCD 2024, Benasque,
Spain, 14-20 January, 2024

J. A. Oller
Departamento de Física, Universidad de Murcia, E-30071 Murcia, Spain

###### Abstract

When a resonance lies near the threshold of a heavier channel, an interesting
feature can occur. The paradigmatic example employed here is the scalar
isoscalar $f_{0}(980)$ resonance that couples to the lighter $\pi\pi$ and
heavier $K\bar{K}$ channels. It is shown that the decay width is given by the
sum or subtraction of the partial decay widths depending on whether the pole
lies in the Riemann sheet that is contiguous with the physical one above or
below the $K\bar{K}$ threshold, respectively. Next, we show that the usually
disregarded renormalization of bare parameters in Flatté or energy-dependent
Breit-Wigner parameterizations is essential to extract physical information.
The compositeness of the $f_{0}(980)$ by using a Flatté parameterization
matched to reproduce the pole properties obtained from Roy equations and other
analytic constraints is evaluated.

##  1 Introduction

A Flatté parameterization [1] is typically used for describing a resonance
that lies near a heavier threshold. Let us denote by $i=1$ and 2 the light and
heavy channels, respectively. To fix ideas think of the $f_{0}(980)$ and the
scalar isoscalar channels $\pi\pi$ and $K\bar{K}$, in this order. Then, around
the $K\bar{K}$ threshold an $S$-wave amplitude is written as

$\displaystyle t_{ij}$ $\displaystyle=$
$\displaystyle\frac{\tilde{g}_{i}\tilde{g}_{j}}{E-E_{f}+i\frac{\widetilde{\Gamma}_{1}}{2}+\frac{i}{2}\tilde{g}_{2}^{2}\sqrt{m_{2}E}}~{},$
(1)

with $E$ the total energy measured with respect to the two-kaon threshold. The
kinematics for the $K\bar{K}$ channel is treated nonrelativistically. This
parameterization is determined by three bare parameters: The bare coupling
$\tilde{g}_{2}$, the bare width $\widetilde{\Gamma}_{1}$ and the bare
resonance mass $E_{f}$. $\widetilde{\Gamma}_{1}$ is related to the bare
coupling $\tilde{g}_{1}$ by

$\displaystyle\widetilde{\Gamma}_{1}$ $\displaystyle=$
$\displaystyle\tilde{g}_{1}^{2}p_{1}~{},$ (2)

where $p_{1}$ is the $\pi\pi$ momentum at the resonance mass $M_{R}=\Re
E_{R}$, and $E_{R}$ is the resonance pole position in Eq. (1). This
parameterization was extensively used in Ref. [2] to study the compositeness
of the $f_{0}(980)$ and $a_{0}(980)$ resonances [3, 4], and revisited in Ref.
[5]. An issue in the traditional way of analyzing Flatté parameterizations was
emphasized in the latter reference. To illustrate it, we gather in Table 1 the
set of Flatté parameterizations for the $f_{0}(980)$ used in [2]. From left to
right, we give the original reference, binding momentum $p_{2}$,
$M_{R}=\Re(p_{2}^{2})/m_{K}$, the pole width $\Gamma=-2\Im p_{2}^{2}/m_{K}$,
and the bare parameters. It is striking that for all cases
$\Gamma\ll\widetilde{\Gamma}_{1}$. Related to that (as shown below), let us
note that all the poles have $\Im p_{2}>0$.

Ref. |  $p_{2}$ (MeV) |  $M_{R}$ (MeV) |  $\Gamma$ (MeV) |  $\widetilde{\Gamma}_{1}$ (MeV) | $\tilde{g}_{2}^{2}$ |  $E_{f}$ (MeV)
---|---|---|---|---|---|---
[6] | $-65+i\,97$ | $981$ | $50.8$ | $149$ | $1.51$ | $-84.3$
[7] | $-58+i\,107$ | $975$ | $50.1$ | $196$ | $2.51$ | $-151.5$
[8] | $-84+i\,17$ | $1005$ | $11.6$ | $129$ | $1.31$ | $+4.6$
[9] | $-69+i\,100$ | $981$ | $55.6$ | $253$ | $2.84$ | $-154$
Table 1: Set of Flatté parameterizations for the $f_{0}(980)$ considered in
[2, 5]. See the text for more details.

Another issue stressed in Ref. [5] concerns the pole determination based on
Roy equations and other analytical constraints in Ref. [10]

$\displaystyle 2m_{K}+E_{R}$ $\displaystyle=$ $\displaystyle 996\pm
7-i\,25^{+10}_{-6}\,\text{MeV}~{},~{}~{}g_{1}=0.46\pm 0.04~{},$ (3)

where $g_{1}$ is the physical coupling to $\pi\pi$ (adopted to our
normalization), obtained from the residue of the partial-wave amplitude at the
$f_{0}(980)$ resonance pole. Let us notice that the partial-decay width to
$\pi\pi$ can be straightforwardly calculated by taking $g_{1}$ from Eq. (3)
into Eq. (2),

$\displaystyle\Gamma_{\pi\pi}=g_{1}^{2}p_{1}=100^{+20}_{-17}\,\text{MeV}~{}\cite[cite]{[\@@bibref{}{Burkert:2022bqo}{}{}]}~{},$
(4)

which is around a factor 2 larger than the pole width $\Gamma$ from Eq. (3).
How can it be?

##  2 Interplay with Riemann sheets

The $f_{0}(980)$ lies close to the two-kaon threshold near 1 GeV.
Consequently, its physical imprint is largely dependent on the Riemann sheet
(RS) in which it lies. We characterize the different RS’s [12] by the signs of
the imaginary parts of the momenta collected as $(\pm,\pm)$, with the 1st(2nd)
sign for $p_{1}(p_{2})$. In this way, $(+,+)$ is the physical or 1st RS,
$(-,+)$ the 2nd RS, $(-,-)$ the third one, and $(+,-)$ the fourth RS. This is
because the square root in the calculation of the momentum as a function of
energy has a right-hand cut. Therefore, when crossing the real energy axis in
between the $\pi\pi$ and $K\bar{K}$ thresholds the 2nd RS $(-,+)$ connects
smoothly with the physical one; when the real axis is crossed above the
$K\bar{K}$ threshold the 3rd RS $(-,-)$ is the one that connects smoothly with
the physical RS.

It is apparent from Table 1 that all the $f_{0}(980)$ poles there are in the
2nd RS, because of the positive sign of $\Im p_{2}$. The pole in Eq. (3) from
Ref. [10] is also in the 2nd RS. This has its importance because of the term
$ig_{2}^{2}p_{2}/2$ in the denominator in Eq. (1). One can consider only one
complex $E$ plane by using the common convention in numerical calculations,
like in Fortran, such that the square root has a left-hand cut. In this way,
the 1st RS corresponds to $\Im E>0$, and Eq. (1) applies. For $\Im E<0$ the
expressions for the different RS’s are:

2nd RS: $\displaystyle t_{ij}=$
$\displaystyle\frac{\tilde{g}_{i}\tilde{g}_{j}}{E-E_{f}+i\frac{\widetilde{\Gamma}_{1}}{2}{\color[rgb]{1,0,0}\mathbf{-}}\frac{i}{2}\tilde{g}_{2}^{2}\sqrt{m_{2}E}}~{},$
(5) 3rd RS: $\displaystyle t_{ij}=$
$\displaystyle\frac{\tilde{g}_{i}\tilde{g}_{j}}{E-E_{f}+i\frac{\widetilde{\Gamma}_{1}}{2}{\color[rgb]{1,0,0}\mathbf{+}}\frac{i}{2}\tilde{g}_{2}^{2}\sqrt{m_{2}E}}~{}.$

We have stressed in Eq. (5) the sign in front of $p_{2}=\sqrt{m_{2}E}$. Then,
we see that in the 3rd RS the $K\bar{K}$ and $\pi\pi$ partial-decay widths add
up, whereas in the 2nd RS $p_{2}$ flips its sign and the $K\bar{K}$
contribution is subtracted to the $\pi\pi$ one to get the pole width. As a
result, we have the following relations between the pole width $\Gamma$ and
the partial-decay widths:

$\displaystyle\text{2nd RS:~{}}\Gamma$
$\displaystyle=\Gamma_{\pi\pi}-\Gamma_{K\bar{K}}~{},$ (6)
$\displaystyle\text{3rd RS:~{}}\Gamma$
$\displaystyle=\Gamma_{\pi\pi}+\Gamma_{K\bar{K}}~{}.$ (7)

The subtraction between the partial-decay widths in Eq. (6) for the pole in
the 2nd RS was first unveiled in Ref. [5]. Two corollaries follow from Eq.
(6): 1) Since $\Gamma>0$ then $\Gamma_{\pi\pi}>\Gamma_{K\bar{K}}$. 2) As
$\Gamma_{\pi\pi}=\Gamma+\Gamma_{K\bar{K}}$ then
$\Gamma_{\pi\pi}>\Gamma_{K\bar{K}}$.

The result in Eq. (6) was then applied in Ref. [11] to the pole position of
the $f_{0}(980)$ in the 2nd RS [10] given in Eq. (3). The fact that
$\Gamma_{\pi\pi}>\Gamma$, cf. Eq. (4), is now understood as due to the
negative contribution of $\Gamma_{K\bar{K}}$ to $\Gamma$. Thus,
$\Gamma_{K\bar{K}}=\Gamma_{\pi\pi}-\Gamma=50^{+26}_{-21}~{}\text{MeV},$ as
calculated in Ref. [11]. Another quantity of interest that was addressed in
this reference is the definition of the total width $\Gamma_{\rm
tot}=\Gamma_{\pi\pi}+\Gamma_{K\bar{K}}$, which does not coincide with the pole
width $\Gamma=\Gamma_{\pi\pi}-\Gamma_{K\bar{K}}$. The definition of
$\Gamma_{\rm tot}$, the same independently of the sheet in which the pole
lies, reflects the fact that in an event distribution the total sum of
resonant events is the sum of events for every channel separately. We report
from Ref. [11] the following results

$\displaystyle\Gamma_{\rm tot}$ $\displaystyle=$
$\displaystyle\Gamma_{\pi\pi}+\Gamma_{K\bar{K}}=151^{+44}_{-37}~{}\text{MeV}~{},$
(8) $\displaystyle\text{BR}_{\pi\pi}$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{\pi\pi}}{\Gamma_{\rm tot}}=0.67\pm 0.07~{},$
$\displaystyle\text{BR}_{K\bar{K}}$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{K\bar{K}}}{\Gamma_{\rm tot}}=0.33\pm 0.07~{},$
$\displaystyle r_{K\bar{K}/\pi\pi}$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{K\bar{K}}}{\Gamma_{\pi\pi}}=0.49\pm 0.11~{}.$

##  3 Renormalization of bare parameters in phenomenological parameterizations

Let us illustrate the process for the Flatté parameterization in Eq. (1),
though the same line of argumentation could be applied to an energy dependent
Breit-Wigner as well. We refer to [5] for the more detailed and original
derivation.

The pole position $E_{R}$ in Eq. (1) near to the $K\bar{K}$ threshold is [5]

$\displaystyle E_{R}$
$\displaystyle=E_{f}-\frac{m_{K}}{8}g_{2}^{4}-\frac{i}{2}\widetilde{\Gamma}_{\pi\pi}+\sigma\frac{g_{2}^{2}}{2}\sqrt{m_{K}\left(\frac{m_{K}g_{2}^{4}}{16}-E_{f}+\frac{i}{2}\widetilde{\Gamma}_{\pi\pi}\right)}~{},$
(9)

where $\sigma=+1(-1)$ corresponds to the pole lying the 2nd(3rd) RS. Next, it
is important to calculate the behavior of the denominator of $t_{ij}(E)$ in
Eq. (1) for $E\to E_{R}$. One has that [5]

$\displaystyle\beta$ $\displaystyle\equiv$ $\displaystyle\lim_{E\to
E_{R}}\frac{E-E_{R}}{E-E_{f}+i\frac{\widetilde{\Gamma}_{1}}{2}+\frac{i}{2}g_{2}^{2}\sqrt{m_{2}E}}$
(10) $\displaystyle=$ $\displaystyle
4\sqrt{|E_{R}|}{m_{k}g_{2}^{4}+16|E_{R}|+4\sigma
g_{2}^{2}\sqrt{2m_{K}(|E_{R}|-M_{R})}}~{},$

with $|E_{R}|=\sqrt{M_{R}^{2}+\Gamma^{2}/4}$. Therefore, the renormalized or
physical couplings in a Flatté parameters are

$\displaystyle g_{i}=\beta^{\frac{1}{2}}\tilde{g}_{i}~{},$ (11)

such that the physical width to $\pi\pi$ is
$\Gamma_{\pi\pi}=\beta\widetilde{\Gamma}_{\pi\pi}$. The renormalized couplings
$g_{i}$ are the ones that must be compared with the couplings obtained by
evaluating the residue of a $T$-matrix, like $g_{1}$ given in Eq. (3) from
Ref. [10]. It is important to stress this point because it is common in the
literature to use the bare couplings and widths of a Flatté parameterization
as physical ones. To see the dramatic impact of $\beta$ we show in Table 2 the
values of $\beta$, $\Gamma_{\pi\pi}=\beta\widetilde{\Gamma}_{\pi\pi}$ and
$\Gamma_{K\bar{K}}=\Gamma_{\pi\pi}-\Gamma$ that correspond to the same Flatté
analyses as in Table 1. It is obvious from the new results in Table 2 that
$\widetilde{\Gamma}_{\pi\pi}$ is very different to the physical one
$\Gamma_{\pi\pi}$, and that $\Gamma_{\pi\pi}>\Gamma$. Thus, Eqs. (6) and (11)
allow to properly extract the physical couplings and widths from a Flatté
parameterization.

Ref. |  $p_{2}$ (MeV) |  $\Gamma$ (MeV) | $\beta$ |  $\Gamma_{\pi\pi}=\beta\widetilde{\Gamma}_{\pi\pi}$ (MeV) | $\Gamma_{K\bar{K}}=\Gamma_{\pi\pi}-\Gamma$
---|---|---|---|---|---
[6] | $-65+i\,97$ | $50.8$ | 0.40 | 59.7 | 8.9
[7] | $-58+i\,107$ | $50.1$ | 0.29 | 56.5 | 6.4
[8] | $-84+i\,17$ | $11.6$ | 0.27 | 67.2 | 11.6
[9] | $-69+i\,100$ | $55.6$ | 0.43 | 55.7 | 44.1
Table 2: The $\beta$ parameter and partial decay widths for the Flatté
parameterizations in Table 1. Notice the smallness of $\beta$.

###  3.1 Compositeness analysis

Now, let us apply the compositeness relation from Refs.[13, 14] to the
$f_{0}(980)$ pole from Ref. [10]

$\displaystyle X$
$\displaystyle=X_{1}+X_{2}~{},~{}X_{1}=\gamma_{1}^{2}\left|\frac{\partial
G_{1}}{\partial
s}\right|_{s_{R}}~{},~{}X_{2}=\gamma_{2}^{2}\left|\frac{\partial
G_{2}}{\partial s}\right|_{s_{R}}~{}.$ (12)

In this equation $X_{i}$ is the partial compositeness of channel $i$ and $X$
is the total compositeness. Let us recall that the a partial compositeness is
the weight of this channel in the composition of the state, and the total
compositeness is the total weight of the meson-meson components. Regarding the
different ingredients in Eq. (12): 1) $s$ is the Mandelstam variable
$s=P^{2}$, with $P$ the total four-momentum, and $s_{R}=(2m_{K}+E_{R})^{2}$.
2) The couplings $\gamma_{i}$ are just proportional to $g_{i}$, such that
$\gamma_{i}=g_{i}\sqrt{8\pi\Re s_{R}}$. 3) The functions $G_{i}(s)$ are the
relativistic unitarity loop functions

$\displaystyle G_{i}$
$\displaystyle=-\frac{1}{16\pi^{2}}\ln\frac{\sigma(s)-1}{\sigma(s)+1}~{},~{}\sigma(s)=\sqrt{1-\frac{4m_{i}^{2}}{s}}~{}.$
(13)

To establish a Flatté parameterization requires three parameters
$(\tilde{g}_{2},\widetilde{\Gamma}_{\pi\pi},E_{f})$. From the pole position
[10] in Eq. (3) we can fix two parameters, but one more is still necessary.
Then, as in Ref. [5], we take $X$ as the third input, and calculate as a
function of it the physical quantities. It turns out that only for $X>0.6$ the
value of $g_{1}$ that results is compatible with Eq. (3). In more detail, we
have $(X,g_{1})=(1,0.47)$, $(0.8,0.45)$, and $(0.6,0.42)$. In the same order,
$X_{2}=0.96$, 0.76 and 0.57, respectively, with $X_{1}\lesssim 0.04\ll X_{2}$
for all cases. Therefore, the application of the compositeness relation of Eq.
12 to a Flatté parameterization required to reproduce the pole properties of
the $f_{0}(980)$ from Ref. [10], drives to the conclusion that this pole is
mainly of composite nature, with its composition dominated by the $K\bar{K}$
contribution.

In summary, we have shown that the pole width for a pole lying in the 2nd RS
is given by the subtraction of the partial-decay widths. We have also
discussed the renormalization process of the bare parameters in Flatté or
energy-dependent Breit-Wigner parameterizations. For clarification and to
avoid confusion, phenomenological analyses should provide the renormalized
parameters, the physically meaningful ones, which can be worked out
straightforwardly from the bare ones directly employed in these
parameterizations.

Acknowledgements. I would like to thank partial financial support to the Grant
PID2022-136510NB-C32 funded by
MCIN/AEI/10.13039/501100011033/ and FEDER, UE.

## References

  * [1] Stanley M. Flatte.  Phys. Lett. B, 63:224–227, 1976.
  * [2] V. Baru, J. Haidenbauer, C. Hanhart, Yu. Kalashnikova, and Alexander Evgenyevich Kudryavtsev.  Phys. Lett. B, 586:53–61, 2004.
  * [3] J. A. Oller, E. Oset, and J. R. Pelaez.  Phys. Rev. D, 62:114017, 2000.
  * [4] M. Albaladejo, J. A. Oller, and L. Roca.  Phys. Rev. D, 82:094019, 2010.
  * [5] Ze-Qiang Wang, Xian-Wei Kang, J. A. Oller, and Lu Zhang.  Phys. Rev. D, 105(7):074016, 2022.
  * [6] R. R. Akhmetshin et al.  Phys. Lett. B, 462:380, 1999.
  * [7] M. N. Achasov et al.  Phys. Lett. B, 485:349–356, 2000.
  * [8] N. N. Achasov and V. V. Gubin.  Phys. Rev. D, 63:094007, 2001.
  * [9] Antonella Antonelli.  eConf, C020620:THAT06, 2002.
  * [10] R. Garcia-Martin, R. Kaminski, J. R. Pelaez, and J. Ruiz de Elvira.  Phys. Rev. Lett., 107:072001, 2011.
  * [11] V. Burkert et al.  Phys. Lett. B, 844:138070, 2023.
  * [12] José Antonio Oller.  Springer Briefs in Physics. Springer, 2019.
  * [13] Zhi-Hui Guo and J. A. Oller.  Phys. Rev. D, 93(9):096001, 2016.
  * [14] J. A. Oller.  Annals Phys., 396:429–458, 2018.