# Random Effect Restricted Mean Survival Time Model

Keisuke Hanada  Kyowa Kirin Co., Ltd  Masahiro Kojima111Address: Biometrics
Department, R&D Division, Kyowa Kirin Co., Ltd. Otemachi Financial City Grand
Cube, 1-9-2 Otemachi, Chiyoda-ku, Tokyo 100-004, Japan. Tel: +81-3-5205-7200
E-Mail<EMAIL_ADDRESS>Kyowa Kirin Co., Ltd

###### Abstract

The restricted mean survival time (RMST) model has been garnering attention as
a way to provide a clinically intuitive measure: the mean survival time. RMST
models, which use methods based on pseudo time-to-event values and inverse
probability censoring weighting, can adjust covariates. However, no approach
has yet been introduced that considers random effects for clusters. In this
paper, we propose a new random-effect RMST. We present two methods of analysis
that consider variable effects by i) using a generalized mixed model with
pseudo-values and ii) integrating the estimated results from the inverse
probability censoring weighting estimating equations for each cluster. We
evaluate our proposed methods through computer simulations. In addition, we
analyze the effect of a mother’s age at birth on under-five deaths in India
using states as clusters.

Keywords: restricted mean survival time, random effect, cluster effect,
pseudo-value, inverse probability censoring weight.

##  1 Introduction

A restricted mean survival time (RMST) is statistically robust and independent
under the assumption of the proportional hazards property, and it gives
clinically interpretable results in terms of mean survival time [1]. When
comparing mean survival times between groups, covariates must be adjusted to
account for bias in the baseline information. To analyze time-to-event data,
how to handle right-censored data should be considered. Unlike Cox regression,
RMST cannot handle a partial likelihood over uncensored data. Andersen et al.
[2] proposed a covariate adjustment method that imputes the pseudo-value (PV)
derived from the leave-one-out method for censored data and uses the
generalized linear model to analyze completed data. A feature of the PV method
is that it can also use censored data. Tian et al. [3] proposed a method using
inverse probability censoring weighting (IPCW) to account for censored data.
There is no approach that considers random effects for clusters.

Considering a random effect in the model allows for different effect sizes
among the clusters, while enabling a common overall effect size. In practice,
treatment effects may differ across clusters, for example, by country or
populations classified by certain biomarker values. Analysis can be performed
to account for heterogeneity across clusters. By estimating the parameters of
the distribution that the variant effect follows, we can also check how much
heterogeneity there may be. In addition, the amount of effect within each
cluster can be accurately estimated according to the degree of heterogeneity
using reduced empirical Bayesian estimation.

In this paper, we propose the random-effect RMST model. We present two methods
of analysis that take into account variable effects by i) using a generalized
mixed model with PVs and ii) integrating the estimated results from the IPCW
estimating equations for each cluster. We evaluate our proposed methods
through computer simulations. In addition, using states as clusters, we
analyze the effect of a mother’s age at birth on under-five deaths in India.

This paper is organized as follows. Section 2 introduces two already proposed
RMST models, proposes the new random-effect RMST model, and introduces the
method used to estimate the model parameter and variance parameter of the
random effect. Section 3 describes the setting and results of the computer
simulations. Section 4 presents the results of an analysis of data from the
eight Empowered Action Group (EAG) states in India. Finally, Section 5
discusses the results.

##  2 Method

We assume that the survival time is $T$, the right-censoring variable is $C$,
and the $p$ covariates are
${\text{\boldmath$X$}}=(X_{1},X_{2},\ldots,X_{p})^{T}$. There are $I$
clusters, the number of subjects in the $i$-th cluster is $N_{i}$, and the
total sample size is $N(=\sum^{I}_{i=1}N_{i})$. The restricted time for
calculating the RMST is a time point $\tau$ of clinical interest. The
observable data of the $n_{ij}$-th subject in cluster $i$ are
$Y_{ij}=\min(T_{ij},C_{ij},\tau)$,
$\Delta_{ij}=I(T_{ij}\leq\min(C_{ij},\tau))$, and
${\text{\boldmath$X$}}_{ij}$.

We propose two types of RMST model with random effects. The first one is based
on a method using PVs proposed by Andersen et al. [4], and we call it the PV
method. The other is based on the IPCW estimating equation proposed by Tian et
al [3] and two-stage estimation method, and we call this method the IPCW
method. We present the first method using PVs. For censored data, PVs are
generated by the leave-one-out method. The RMST estimator for all subjects is
${\hat{m}}(\tau)=\int^{\tau}_{0}{\hat{S}}(t)dt$, where ${\hat{S}}(t)$ is the
Kaplan-Meier estimator. The RMST estimator excluding subject $n_{ij}$ is
${\hat{m}}_{-ij}(\tau)$, and the PV of subject $n_{ij}$ is
$N{\hat{m}}(\tau)-(N-1){\hat{m}}_{-ij}(\tau)$. In this paper, we compute the
PVs within each cluster to account for the heterogeneity among clusters. The
RMST estimator in the $i$-th cluster is

$\displaystyle{\hat{m}}^{(i)}(\tau)=\int^{\tau}_{0}{\hat{S}}^{(i)}(t)dt,$ (1)

and the PVs are calculated from

$\displaystyle
N_{i}{\hat{m}}^{(i)}(\tau)-(N_{i}-1){\hat{m}}^{(i)}_{-ij}(\tau).$ (2)

Because all censored data are complemented by PVs, the RMST can be calculated
by a model analysis without considering the censoring. We consider an analysis
that accounts for a random effect using a generalized linear mixed model with
the link function as
$g(\mu_{ij})={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$x$}}_{ij}+v_{i}$,
where $\mu_{ij}=E[Y_{ij}|v_{i}]$ and $v_{i}$ is a random effect distributed as
a normal distribution of ${\cal N}(0,{\sigma}^{2}_{v})$. The parameters of
$\beta$ and ${\sigma}^{2}_{v}$ can be estimated by a generalized estimating
equation, such as Laplace approximation, adaptive Gaussian-Hermite quadrature,
or penalized quasi-likelihood. Note that we introduce the generalized
estimating equation method by referring to Jiang and Nguyen[5]. We denote the
observed time-to-event values and the generated PVs as $Y$ without
distinction. We assume the expectation of $Y$ is

$\displaystyle E[Y_{ij}]=\int
h({\text{\boldmath$x$}}_{ij}^{T}{\text{\boldmath$\beta$}}+{\sigma}_{v}u_{i})f(u_{i})du_{i}\equiv
M_{i}({\text{\boldmath$\theta$}}),$ (3)

where $h(\cdot)=g^{-1}(\cdot)$, $u_{i}$ is a standard normal variable,
$f(u_{i})$ is a standard normal density function, and
${\text{\boldmath$\theta$}}=({\text{\boldmath$\beta$}}^{T},{\sigma}_{v})^{T}$.
The first derivatives are

$\displaystyle\frac{\partial
M_{i}({\text{\boldmath$\theta$}})}{\partial{\text{\boldmath$\theta$}}}=\begin{pmatrix}\frac{\partial
M_{i}({\text{\boldmath$\theta$}})}{\partial{\text{\boldmath$\beta$}}}\\\
\frac{\partial M_{i}({\text{\boldmath$\theta$}})}{\partial{\sigma}_{v}}\\\
\end{pmatrix}=\begin{pmatrix}{\text{\boldmath$x$}}_{ij}\int
h^{\prime}({\text{\boldmath$x$}}_{ij}^{T}{\text{\boldmath$\beta$}}+{\sigma}_{v}u_{i})f(u_{i})du_{i}\\\
\int
h^{\prime}({\text{\boldmath$x$}}_{ij}^{T}{\text{\boldmath$\beta$}}+{\sigma}_{v}u_{i})u_{i}f(u_{i})du_{i}\\\
\end{pmatrix}.$ (4)

We assume $V[Y_{ij}]\equiv V_{0}$. The generalized estimating equation is

$\displaystyle\sum^{I}_{i=1}\frac{\partial{\text{\boldmath$M$}}_{i}({\text{\boldmath$\theta$}})^{T}}{\partial{\text{\boldmath$\theta$}}}{\widehat{V}}_{i}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$M$}}_{i}({\text{\boldmath$\theta$}}))=0,$
(5)

where ${\text{\boldmath$V$}}_{i}$ is a working covariance matrix for the
$i$-th cluster. The initial value of ${\text{\boldmath$V$}}_{i}$ is an
$N_{i}\times N_{i}$ identity matrix ${\text{\boldmath$V$}}_{i}$. Parameters
$\theta$ are estimated from Equation (4). Next, we calculate
${\text{\boldmath$V$}}_{i}$ using the estimated values again and iterate
through Equations (4) and (5). If there is no difference between the current
estimation values and the previous estimation values, the iterations are
stopped, and the last estimation values are output as the generalized
estimating equation estimators.
From here, we introduce the shrinkage estimator and adjusted confidence
interval for each cluster using the empirical Bayes estimator when $Y$ is
distributed as a normal distribution [6]. Specifically, we are interested in
the difference between groups. We assume the indicator variable for the group
is $x_{1ij}$ and the parameter of the group is ${\beta}_{1}$. The empirical
Bayes estimator of ${\beta}_{1i}$ for the $i$-th cluster is

$\displaystyle{\hat{\beta}}_{1i}^{EB}={\hat{\beta}}_{1i}+\frac{{\hat{\sigma}}_{i}^{2}}{{\hat{\sigma}}_{1i}^{2}+{\hat{\sigma}}_{v}^{2}}({\hat{\beta}}_{1}-{\hat{\beta}}_{1i}),$
(6)

where ${\hat{\sigma}}^{2}$ is the square of the standard error of
${\hat{\beta}}_{1}$,
${\overline{x}}_{1i}=\frac{1}{N_{i}}\sum^{N_{i}}_{j=1}x_{1ij}$, and
${\overline{y}}_{1i}=\frac{1}{N_{i}}\sum^{N_{i}}_{j=1}y_{1ij}$. The adjusted
$100(1-{\alpha})\%$ confidence interval is
$\left[{\hat{\beta}}_{1i}^{EB}-z_{{\alpha}/2}\sqrt{\frac{{\hat{\sigma}}^{2}{\hat{\sigma}}_{v}^{2}}{{\hat{\sigma}}^{2}+{\hat{\sigma}}_{v}^{2}}},{\hat{\beta}}_{1i}^{EB}+z_{{\alpha}/2}\sqrt{\frac{{\hat{\sigma}}^{2}{\hat{\sigma}}_{v}^{2}}{{\hat{\sigma}}^{2}+{\hat{\sigma}}_{v}^{2}}}\right]$.
If $Y$ is distributed as a log-normal distribution, then we refer the reader
to Berg and Chandra [7]. For the log-normal distribution, the formula for the
error in the confidence interval is very complicated.

The second proposed method is based on the IPCW estimating equation by Tian et
al. [3]. Because the IPCW estimating equation is unable to account for a
random effect, we consider a two-stage estimation, first estimating the
coefficient parameters by cluster and then accounting for heterogeneity across
clusters in the estimated results. For the $i$-th cluster, the parameters are
estimated by the following IPCW estimating equation:

$\displaystyle\frac{1}{N_{i}}\sum^{N_{i}}_{j=1}\frac{{\widehat{\Delta}}_{ij}}{KM(y_{ij})}{\text{\boldmath$x$}}_{ij}(y_{ij}-h({\text{\boldmath$x$}}_{ij}^{T}{\text{\boldmath$\beta$}}))=0,$
(8)

where ${\widehat{\Delta}}_{ij}=I(Y_{ij}\leq C_{ij})$ and $KM(y_{ij})$ is the
Kaplan-Meier estimator of the censoring time $C_{ij}$. We assume that the
estimator of the parameters in the $i$-th cluster is
${\widehat{\text{\boldmath$\beta$}}}_{i}$. We assume that the variable of our
particular interest is the $1$-th parameter of the group, and the standard
error of ${\hat{\beta}}_{1i}$ is ${\hat{\sigma}}_{1i}$. There is a random
effect $v\sim{\cal N}(0,{\sigma}^{2}_{v})$. The estimator of
${\sigma}^{2}_{v}$ is

$\displaystyle{\hat{\sigma}}^{2}_{v}=\frac{Q-(I-1)}{c},$ (9)

where
$Q=\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-2}({\hat{\beta}}_{1i}-{\overline{\hat{\beta}}}_{1})^{2}$,

$\displaystyle{\overline{\hat{\beta}}}_{1}=\frac{\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-2}{\hat{\beta}}_{1i}}{\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-2}}\mbox{,
and
}c=\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-2}-\frac{\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-4}}{\sum^{I}_{i=1}{\hat{\sigma}}_{1i}^{-2}}.$
(10)

Here, ${\overline{\hat{\beta}}}_{1}$ is the pooled mean estimator without a
random effect. If $Q<I-1$, then ${\hat{\sigma}}^{2}_{v}=0$. The pooled mean
estimator ${\overline{\overline{\hat{\beta}}}}_{1}$ with a random effect is
$\sum^{I}_{i=1}({\hat{\sigma}}_{1i}^{2}+{\hat{\sigma}}_{v}^{2})^{-1}{\hat{\beta}}_{1i}/\sum^{I}_{i=1}({\hat{\sigma}}_{1i}^{2}+{\hat{\sigma}}_{v}^{2})^{-1}$.
Each shrinkage estimator is
${\tilde{\beta}}^{EB}_{1i}={\hat{\beta}}_{1i}+\frac{{\hat{\sigma}}_{i}^{2}}{{\hat{\sigma}}_{1i}^{2}+{\hat{\sigma}}_{v}^{2}}({\overline{\overline{\hat{\beta}}}}_{1}-{\hat{\beta}}_{1i})$
[8]. The adjusted $100(1-{\alpha})\%$ confidence interval is
$\left[{\tilde{\beta}}^{EB}_{1i}-z_{{\alpha}/2}\sqrt{\frac{{\hat{\sigma}}^{2}_{v}{\hat{\sigma}}^{2}_{1i}}{{\hat{\sigma}}^{2}_{v}+{\hat{\sigma}}^{2}_{1i}}},{\tilde{\beta}}^{EB}_{1i}+z_{{\alpha}/2}\sqrt{\frac{{\hat{\sigma}}^{2}_{v}{\hat{\sigma}}^{2}_{1i}}{{\hat{\sigma}}^{2}_{v}+{\hat{\sigma}}^{2}_{1i}}}\right]$
[9]. For the log-normal distribution, we refer the reader to Slud and Maiti
(2006) [10].

##  3 Simulation

We evaluated the performance of our proposal method using Monte Carlo
simulation. Data-generating model 1 is

$\displaystyle
Y_{ij}=\exp({\beta}_{0}+{\beta}_{1}x_{1ij}+v_{i}){\varepsilon}_{ij},$ (11)

where ${\beta}_{0}=1$, ${\beta}_{1}=0.5$, $v_{i}\sim{\cal
N}(0,{\sigma}^{2}_{v})$, ${\sigma}^{2}_{v}=0.3^{2}$, and
${\varepsilon}_{ij}\sim Exp(1)$. Moreover, $x_{1ij}$ is 1 for the treatment
group and 0 for the control group. The groups are assigned randomly and
uniformly. In addition, we consider a model with covariates. Data-generating
model 2 is

$\displaystyle
Y_{ij}=\exp({\beta}_{0}+{\beta}_{1}x_{1ij}+{\beta}_{2}x_{2ij}+{\beta}_{3}x_{3ij}+v_{i}){\varepsilon}_{ij},$
(12)

where ${\beta}_{0}=1$, ${\beta}_{1}=0.5$, ${\beta}_{2}=0.1$,
${\beta}_{3}=-0.5$, $v_{i}\sim{\cal N}(0,{\sigma}^{2}_{v})$,
${\sigma}^{2}_{v}=0.3^{2}$, and ${\varepsilon}_{ij}\sim Exp(1)$. Furthermore,
$x_{1ij}$ is 1 for the treatment group and 0 for the control group, $x_{2ij}$
is ${\cal N}(1,0.5^{2})$, and $x_{3ij}$ is set to 1 or 0 randomly. The groups
are assigned randomly and uniformly. The censoring probability is set to
$0.1$, $0.5$, and $0.9$. We consider the three numbers of clusters: 5, 8, and
10. The number of simulations is 10,000, and the total sample size is 400 or
1,000. The restricted time $\tau$ is 5 months. The evaluation points are

  1. 1.

Average bias (the estimated result - true value) of $\beta_{1}$;

  2. 2.

Mean squared error (MSE) of $\beta_{1}$;

  3. 3.

Coverage probability of the confidence interval of $\beta_{1}$;

  4. 4.

Length of the confidence interval of $\beta_{1}$;

  5. 5.

Average bias of ${\sigma}^{2}_{v}$;

  6. 6.

MSE of ${\sigma}^{2}_{v}$.

The reference model for comparing performance is a model that excludes the
variable effect $v_{i}$ from Equations (11) and (12). The estimator of
parameter $\hat{\beta}_{1}$ from the generating model may not converge to the
true value of $\beta_{1}$. Therefore, we assume that the true parameter value
for the treatment group of the RMST model is $\hat{\beta}_{1}$. Moreover, we
calculate the approximate true value of parameter $\bar{\beta}_{1}$ using the
100,000 data points for each group. The true value $\bar{\beta}_{1}$ is given
in Appendix A.2. To evaluate the coverage probability of the confidence
interval of $\beta_{1}$, we calculate the $95\%$ confidence interval as
$[{\overline{\hat{\beta}}}_{1}-z_{97.5}\sum_{i=1}^{I}\hat{\sigma}_{i}^{-2},{\overline{\hat{\beta}}}_{1}+z_{97.5}\sum_{i=1}^{I}\hat{\sigma}_{i}^{-2}]$.

The simulation results with censoring probabilities of $0.1$ and $0.5$ are
presented in Tables 1 and 2. The result with a censoring probability of $0.9$
is given in Appendix A.1.

Average bias of $\beta_{1}$: The bias of the PV method is generally small in
all situations, but the bias increases with the number of clusters when the
censoring probability is large (0.9) regardless of the model. The bias of the
IPCW method is also small for model 1 and when the censoring probability is
small (0.1). However, under conditions with moderate to large censoring
probabilities, the bias of the IPCW method increases significantly more than
the bias of the PV method when the number of clusters is increased and the
number of subjects is small.

MSE of $\beta_{1}$: The MSE of the PV method decreases as the number of
subjects is increased for any model and censoring probability. The MSE
decreases under a large censoring probability of 0.9, but the correction due
to PVs may have an impact on it. The MSE of the IPCW method also decreases as
the number of subjects is increased, but the MSE of the IPCW method is higher
than that of the PV method because of the larger censoring probability.

Coverage probability of the confidence interval of $\beta_{1}$: The coverage
probability of the PV method approaches a nominal level as the number of
subjects increases. Increasing the number of clusters improves coverage
probability if the number of subjects is sufficiently large, but decreases it
if the number of subjects is not large enough. The IPCW method shows a similar
trend to the PV method, but the effect of the censoring probability is more
significant. In particular, the coverage probability falls below 0.5 when the
number of clusters is increased with a small number of subjects under a large
censoring probability of 0.9.

Length of the confidence interval of $\beta_{1}$: The length of the confidence
interval tends to be the same for both methods. That is, an increase in the
number of subjects shortens the length, but an increase in the number of
clusters does not change the length much. Under moderate to large censoring
probabilities, the confidence interval of model 2 tends to be longer than that
of model 1.

Average bias of ${\sigma}^{2}_{v}$: The bias of the PV method decreases
slightly as the number of clusters increases, but increases as the number of
subjects increases. The result for the IPCW method is similar to that of the
PV method, with bias decreasing with an increase in the number of clusters and
increasing with an increase in the number of subjects. However, the PV method
has a positive bias, in contrast to the IPCW method, which has a negative
bias. For medium to large censoring probabilities and a small number of
subjects, the bias of the IPCW method increases even if the number of clusters
increases.

MSE of ${\sigma}^{2}_{v}$: The MSE of the PV method decreases as the number of
subjects or clusters are increased. The MSE of the IPCW method remains smaller
than that of the PV method in situations where the number of subjects is large
and the censoring probability is small. However, when the censoring
probability is moderate to large, the MSE of the IPCW method increases as the
number of clusters increases.

##  4 Actual data analysis

We re-analyzed the survival data of children up to age 5 years in India during
the years 2019–2021. In India, eight EAG states were established to facilitate
regional programs in areas that have lagged behind in controlling population
growth to manageable levels. We focused on the survival data of the EAG
states. The total sample size is 101,437, which consists of 37,033 individuals
born to mothers 12–19 years old (12–19 group), 63,380 individuals born to
mothers 20–30 years old (20–30 group), and 1,024 individuals born to mothers
31 years and older (31+ group). The total number of deaths is 4,883,
comprising 1,920, 2,889, and 74 deaths for those born to mothers in the 12–19,
20–30, and 31+ groups, respectively. Moreover, the censoring probabilities are
94.8%, 95.4%, and 92.8%, respectively. The baseline information consists of
the sex, place of birth, and order of birth. The number of deaths and censors
for each state are shown on the left side in Figures 1, 2, 3, and 4. We also
re-analyzed data of 2019–2021 of an Indian national family health survey that
included data for children up to age 5 years. The restricted time, which is
the time at which we can identify more than 80% of the deaths (i.e., 4,094 out
of 4,883 total deaths), is 50 months.

The estimated random effects ${\hat{\sigma}}_{v}^{2}$ of the PV method are
0.026 for 12–19 group vs 20–30 group and 0.036 for 12–19 group vs 31+ group.
The estimated random effects ${\hat{\sigma}}_{v}^{2}$ of IPCW method are 0.056
for 12–19 group vs 20–30 group and 0.414 for 12–19 group vs 31+ group. The
RMST differences are shown in Figures 1, 2, 3, and 4.

For the PV method, when random effects are accounted for, the results
comparing the 12–19 group with the 20–30 group are reduced with respect to the
overall results because of the small variance among the states. In particular,
in Chhattisgarh, where the RMST difference deviates significantly from the
overall result, the reduction was larger because the small sample size was
found to bias the results by chance. For the 12–19 group versus the 31+ group,
because of the small variance among states and the much smaller sample size
for the 31+ group, the results for each state obtained by the random-effect
RMST are substantially reduced. This means that in all states, the prognosis
is better if the child is born to a mother who is 31 years or older.

For the IPCW method, the results comparing the 12–19 group with the 20–30
group reveal that, as for the PV method, the variance among states is small;
hence, results that deviate from the overall results, such as the result for
Chhattisgarh, are reduced with respect to the overall results. In addition,
when the 95% confidence intervals are too long because of the small sample
size, as in the result for Uttarakhand, the 95% confidence intervals are
shorter. For the 12–19 group vs 31+ group, because the IPCW method found a
greater variance among the states than the PV method, the 95% confidence
interval was not shortened as much as it was by the PV method, but the 95%
confidence interval could be shortened for results from states with long
confidence intervals.

##  5 Discussion

We proposed two novel random-effect RMST models. One imputes PVs for censored
data to create a complete dataset and use a generalized linear mixed model.
The other estimates the variance parameter of the random effect among clusters
using the estimation results obtained by the IPCW estimating equation for each
cluster. The real data analysis revealed that the shrinkage effect using the
PV method is stronger than that obtained using the IPCW method because the
number of interest events increases when the PVs of censored data are imputed.

The random-effect RMST model using the PV method analyzes all data by imputing
the PVs for censored data. Given a large sample size, the individual PVs are
approximately independent [4], and the estimators obtained from the
generalized linear mixed model are consistent [11, 12]. However, a limitation
is that the PVs are not imputed while taking into account the random effects.
For the IPCW method, because an estimated equation cannot be constructed that
considers the random effects, we analyzed within each cluster and estimated
the random effects of the clusters. The estimator is consistent with respect
to the true parameter value [9, 10]. The limitation of this approach is that
it is not possible to perform an analysis that adjusts for random effects and
covariates using data from all clusters.

The simulation study showed that as the sample size increased, the MSE
decreased and the estimation accuracy improved. For the PV method, as the
censoring probability increased, the accuracy of the estimated coefficient
parameters did not change, but the results of the estimated parameters of the
variance of the random effects decreased. The accuracy of the estimated random
effect increased because of the smaller observation error for the PVs for
censored data. Moreover, the results remained stable despite the increased
censoring probability because the PVs were working effectively. For the IPCW
method, the accuracy of the coefficient parameter estimates decreased as the
censoring probability increased because the number of observed events was
reduced. In addition, the estimated results for the parameter of variance of
the random effect grew worse as the number of clusters increased. However,
increasing the sample size increased the accuracy of estimation. With both
methods, the coverage probability approached the true value of 0.95 as the
accuracy of the coefficient parameter estimation increased.

From the actual study, for the 12–19 group vs 31+ group comparison, when the
number of cases was small, the variance between states differed in two ways.
We believe that this result reflects the results of the simulated study, since
the IPCW method is less accurate in estimation when the number of clusters is
large and the sample size is small. The PV method, which stores the censoring
as a PV, analyzes a larger number of data as time-to-event, and hence the
confidence intervals were shorter than those of the IPCW method, regardless of
the presence of random effects.

We recommend the random-effect RMST model based on the PV method because it
exhibited less bias in the simulation results. This is because the 95%
confidence interval was closer to the 95% probability of containing the true
value.

## References

  * [1] Hajime Uno, Brian Claggett, Lu Tian, Eisuke Inoue, Paul Gallo, Toshio Miyata, Deborah Schrag, Masahiro Takeuchi, Yoshiaki Uyama, Lihui Zhao, et al.  Moving beyond the hazard ratio in quantifying the between-group difference in survival analysis.  Journal of clinical Oncology, 32(22):2380, 2014.
  * [2] Per Kragh Andersen, Mette Gerster Hansen, and John P Klein.  Regression analysis of restricted mean survival time based on pseudo-observations.  Lifetime data analysis, 10:335–350, 2004.
  * [3] Lu Tian, Lihui Zhao, and LJ Wei.  Predicting the restricted mean event time with the subject’s baseline covariates in survival analysis.  Biostatistics, 15(2):222–233, 2014.
  * [4] Per Kragh Andersen, John P Klein, and Susanne Rosth$o$j.  Generalised linear models for correlated pseudo-observations, with applications to multi-state models.  Biometrika, 90(1):15–27, 2003.
  * [5] Jiming Jiang and Thuan Nguyen.  Linear and generalized linear mixed models and their applications, volume 1.  Springer, 2007.
  * [6] Tatsuya Kubokawa.  Corrected empirical bayes confidence intervals in nested error regression models.  Journal of the Korean Statistical Society, 39(2):221–236, 2010.
  * [7] Emily Berg and Hukum Chandra.  Small area prediction for a unit-level lognormal model.  Computational Statistics & Data Analysis, 78:159–175, 2014.
  * [8] Hui Quan, Mingyu Li, Weichung Joe Shih, Soo Peter Ouyang, Joshua Chen, Ji Zhang, and Peng-Liang Zhao.  Empirical shrinkage estimator for consistency assessment of treatment effects in multi-regional clinical trials.  Statistics in Medicine, 32(10):1691–1706, 2013.
  * [9] Masayo Yoshimori and Partha Lahiri.  A second-order efficient empirical bayes confidence interval.  2014\.
  * [10] Eric V Slud and Tapabrata Maiti.  Mean-squared error estimation in transformed fay–herriot models.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 68(2):239–257, 2006.
  * [11] Jiming Jiang.  Consistent estimators in generalized linear mixed models.  Journal of the American Statistical Association, 93(442):720–729, 1998.
  * [12] Jiming Jiang.  The subset argument and consistency of mle in glmm: Answer to an open problem and beyond.  2013\.

Figure 1: (PV method) Forest plot of the difference between the restricted mean survival times (12–19 group vs 20–30 group) Figure 2: (PV method) Forest plot of the difference between the restricted mean survival times (12–19 group vs 31+ group) Figure 3: (IPCW method) Forest plot of the difference between the restricted mean survival times (12–19 group vs 20–30 group) Figure 4: (IPCW method) Forest plot of the difference between the restricted mean survival times (12–19 group vs 31+ group) Table 1: Simulation results of model 1. | Censoring | PV method | IPCW method
---|---|---|---
Cluster | probability | Total N=400 | Total N=1000 | Total N=400 | Total N=1000
Bias ($\beta_{1}$)
5 | 0.1 | -0.010 | -0.009 | -0.008 | -0.008
8 | 0.1 | -0.002 | -0.008 | 0.004 | -0.005
10 | 0.1 | 0.010 | -0.008 | 0.019 | -0.005
5 | 0.5 | -0.014 | -0.020 | -0.003 | -0.028
8 | 0.5 | -0.012 | -0.020 | 0.019 | -0.024
10 | 0.5 | -0.012 | -0.016 | 0.038 | -0.015
MSE ($\beta_{1}$)
5 | 0.1 | 0.034 | 0.013 | 0.036 | 0.014
8 | 0.1 | 0.034 | 0.013 | 0.036 | 0.014
10 | 0.1 | 0.037 | 0.013 | 0.040 | 0.014
5 | 0.5 | 0.039 | 0.016 | 0.060 | 0.025
8 | 0.5 | 0.039 | 0.016 | 0.064 | 0.024
10 | 0.5 | 0.041 | 0.015 | 0.073 | 0.023
Coverage probability of the confidence interval of $\beta_{1}$
5 | 0.1 | 0.936 | 0.940 | 0.936 | 0.940
8 | 0.1 | 0.936 | 0.939 | 0.928 | 0.938
10 | 0.1 | 0.920 | 0.940 | 0.912 | 0.936
5 | 0.5 | 0.932 | 0.931 | 0.885 | 0.879
8 | 0.5 | 0.929 | 0.934 | 0.871 | 0.891
10 | 0.5 | 0.924 | 0.934 | 0.838 | 0.890
Length of the confidence interval of $\beta_{1}$
5 | 0.1 | 0.685 | 0.436 | 0.698 | 0.444
8 | 0.1 | 0.679 | 0.435 | 0.690 | 0.443
10 | 0.1 | 0.675 | 0.434 | 0.685 | 0.442
5 | 0.5 | 0.725 | 0.461 | 0.777 | 0.498
8 | 0.5 | 0.718 | 0.459 | 0.759 | 0.494
10 | 0.5 | 0.712 | 0.459 | 0.745 | 0.492
Bias ($\sigma_{v}^{2}$)
5 | 0.1 | 0.067 | 0.074 | -0.040 | -0.069
8 | 0.1 | 0.066 | 0.074 | -0.018 | -0.064
10 | 0.1 | 0.064 | 0.072 | -0.002 | -0.060
5 | 0.5 | 0.003 | 0.013 | 0.018 | -0.039
8 | 0.5 | 0.001 | 0.012 | 0.085 | -0.026
10 | 0.5 | 0.001 | 0.010 | 0.155 | -0.017
MSE ($\sigma_{v}^{2}$)
5 | 0.1 | 0.025 | 0.021 | 0.010 | 0.006
8 | 0.1 | 0.019 | 0.016 | 0.014 | 0.006
10 | 0.1 | 0.018 | 0.014 | 0.018 | 0.006
5 | 0.5 | 0.013 | 0.009 | 0.031 | 0.007
8 | 0.5 | 0.011 | 0.006 | 0.067 | 0.008
10 | 0.5 | 0.010 | 0.005 | 0.125 | 0.009
Table 2: Simulation results of model 2. | Censoring | PV method | IPCW method
---|---|---|---
Cluster | probability | Total N=400 | Total N=1000 | Total N=400 | Total N=1000
Bias ($\beta_{1}$)
5 | 0.1 | -0.015 | -0.020 | -0.012 | -0.020
8 | 0.1 | -0.009 | -0.017 | -0.003 | -0.015
10 | 0.1 | -0.003 | -0.016 | 0.004 | -0.013
5 | 0.5 | 0.005 | 0.005 | 0.031 | 0.006
8 | 0.5 | 0.010 | 0.004 | 0.062 | 0.011
10 | 0.5 | 0.018 | 0.006 | 0.097 | 0.021
MSE ($\beta_{1}$)
5 | 0.1 | 0.033 | 0.013 | 0.035 | 0.014
8 | 0.1 | 0.034 | 0.013 | 0.036 | 0.014
10 | 0.1 | 0.035 | 0.013 | 0.038 | 0.014
5 | 0.5 | 0.042 | 0.017 | 0.069 | 0.027
8 | 0.5 | 0.042 | 0.016 | 0.079 | 0.026
10 | 0.5 | 0.044 | 0.016 | 0.095 | 0.026
Coverage probability of the confidence interval of $\beta_{1}$
5 | 0.1 | 0.941 | 0.941 | 0.936 | 0.939
8 | 0.1 | 0.934 | 0.942 | 0.931 | 0.941
10 | 0.1 | 0.924 | 0.942 | 0.919 | 0.942
5 | 0.5 | 0.935 | 0.938 | 0.882 | 0.894
8 | 0.5 | 0.929 | 0.938 | 0.850 | 0.899
10 | 0.5 | 0.922 | 0.941 | 0.802 | 0.897
Length of the confidence interval of $\beta_{1}$
5 | 0.1 | 0.683 | 0.434 | 0.697 | 0.444
8 | 0.1 | 0.677 | 0.433 | 0.690 | 0.443
10 | 0.1 | 0.672 | 0.432 | 0.684 | 0.442
5 | 0.5 | 0.759 | 0.482 | 0.825 | 0.530
8 | 0.5 | 0.753 | 0.481 | 0.804 | 0.526
10 | 0.5 | 0.748 | 0.480 | 0.787 | 0.523
Bias ($\sigma_{v}^{2}$)
5 | 0.1 | 0.065 | 0.072 | -0.038 | -0.071
8 | 0.1 | 0.064 | 0.072 | -0.019 | -0.065
10 | 0.1 | 0.063 | 0.072 | -0.002 | -0.060
5 | 0.5 | 0.012 | 0.021 | 0.032 | -0.040
8 | 0.5 | 0.008 | 0.019 | 0.111 | -0.023
10 | 0.5 | 0.009 | 0.019 | 0.195 | -0.011
MSE ($\sigma_{v}^{2}$)
5 | 0.1 | 0.024 | 0.021 | 0.010 | 0.006
8 | 0.1 | 0.019 | 0.015 | 0.014 | 0.006
10 | 0.1 | 0.017 | 0.013 | 0.018 | 0.006
5 | 0.5 | 0.015 | 0.010 | 0.039 | 0.008
8 | 0.5 | 0.012 | 0.007 | 0.092 | 0.009
10 | 0.5 | 0.012 | 0.006 | 0.167 | 0.011

##  Appendix A Appendix

###  A.1 Simulated results for a censoring probability of 0.9

Table 3: Simulation results for a censoring probability of 0.9.

|  | PV method | IPCW method
---|---|---|---
Cluster | Model | Total N=400 | Total N=1000 | Total N=400 | Total N=1000
Bias ($\beta_{1}$)
5 | 1 | -0.018 | -0.006 | 0.023 | -0.008
8 | 1 | -0.052 | -0.013 | 0.056 | -0.010
10 | 1 | -0.082 | -0.019 | 0.039 | -0.008
5 | 2 | -0.019 | -0.005 | 0.049 | -0.008
8 | 2 | -0.050 | -0.014 | 0.135 | -0.001
10 | 2 | -0.082 | -0.018 | 0.152 | 0.013
MSE ($\beta_{1}$)
5 | 1 | 0.015 | 0.006 | 0.098 | 0.016
8 | 1 | 0.016 | 0.006 | 0.322 | 0.022
10 | 1 | 0.018 | 0.006 | 0.431 | 0.029
5 | 2 | 0.018 | 0.008 | 0.151 | 0.022
8 | 2 | 0.019 | 0.007 | 0.463 | 0.036
10 | 2 | 0.022 | 0.008 | 0.779 | 0.056
Coverage probability of the confidence interval of $\beta_{1}$
5 | 1 | 0.890 | 0.926 | 0.717 | 0.797
8 | 1 | 0.713 | 0.918 | 0.480 | 0.758
10 | 1 | 0.502 | 0.897 | 0.316 | 0.708
5 | 2 | 0.905 | 0.931 | 0.686 | 0.780
8 | 2 | 0.758 | 0.923 | 0.453 | 0.722
10 | 2 | 0.563 | 0.910 | 0.307 | 0.680
Length of the confidence interval of $\beta_{1}$
5 | 1 | 0.423 | 0.288 | 0.437 | 0.304
8 | 1 | 0.330 | 0.278 | 0.330 | 0.290
10 | 1 | 0.230 | 0.270 | 0.241 | 0.280
5 | 2 | 0.474 | 0.319 | 0.495 | 0.342
8 | 2 | 0.385 | 0.309 | 0.383 | 0.327
10 | 2 | 0.283 | 0.301 | 0.291 | 0.316
Bias ($\sigma_{v}^{2}$)
5 | 1 | -0.087 | -0.083 | 0.156 | -0.032
8 | 1 | -0.089 | -0.084 | 0.462 | 0.015
10 | 1 | -0.088 | -0.084 | 0.595 | 0.063
5 | 2 | -0.086 | -0.082 | 0.291 | -0.012
8 | 2 | -0.088 | -0.083 | 0.723 | 0.081
10 | 2 | -0.086 | -0.083 | 0.952 | 0.166
MSE ($\sigma_{v}^{2}$)
5 | 1 | 0.008 | 0.007 | 0.483 | 0.026
8 | 1 | 0.009 | 0.007 | 1.798 | 0.079
10 | 1 | 0.009 | 0.007 | 2.726 | 0.144
5 | 2 | 0.009 | 0.007 | 0.907 | 0.043
8 | 2 | 0.009 | 0.007 | 2.557 | 0.173
10 | 2 | 0.010 | 0.007 | 4.680 | 0.327

###  A.2 Simulated true value of parameter $\beta_{1}$

Table 4: Simulated true value of parameter $\bar{\beta}_{1}$ Model | Method |  | Censoring
---
probability
| Number of
---
simulations
$\bar{\beta}_{1}$
Model 1 | PV | 0.1 | 100000 | 0.7156
|  | 0.5 | 100000 | 0.5681
|  | 0.9 | 100000 | 0.1531
| IPCW | 0.1 | 100000 | 0.7623
|  | 0.5 | 100000 | 0.8086
|  | 0.9 | 100000 | 0.3057
Model 2 | PV | 0.1 | 100000 | 0.7156
|  | 0.5 | 100000 | 0.5681
|  | 0.9 | 100000 | 0.1682
| IPCW | 0.1 | 100000 | 0.7623
|  | 0.5 | 100000 | 0.8086
|  | 0.9 | 100000 | 0.3566


#  BPS deformations of the Skyrme model

J. M. Queiruga<EMAIL_ADDRESS>Department of Theoretical
Physics, UPV/EHU, 48080, Bilbao, Spain

(29-11-2020)

###### Abstract

We study several deformations of the Skyrme model in three dimensions with
self-dual sectors of arbitrary baryonic charge. We show that, for a family of
background metrics as well as for a family of field dependent couplings, the
model has one BPS sector, which may have any topological charge. We also study
the gravitating case, where there are infinite BPS sectors, provided that a
cosmological constant is added to the model.

##  1 Introduction

Although the details still remain to be fully understood, the Skyrme model
Skyrme1 ; Skyrme2  can be interpreted as a low energy limit of QCD when the
number of quark colors tends to infinity witten1 ; witten2 . In particular, it
can be used to describe some properties of nuclei, like for example binding
energies. In this effective approach the fundamental fields are interpreted as
pions, while the topological solutions of the model, called Skyrmions,
represent nuclei. From a mathematical point of view, a Skyrmion is a map from
a three dimensional manifold (a one point compactification of $\mathbb{R}^{3}$
in the standard formulation) to $SU(2)$. Since $\pi_{3}(SU(2))=\mathbb{Z}$ the
solutions are classified by a topological number, normally identified with the
baryonic charge, and therefore they are topologically stable.

Despite of its rich mathematical structure, the original formulation of the
Skyrme model only provides a rough description of some properties of the low
energy QCD Nappi . For example, regarding the applications to nuclear physics,
one of the main shortcomings of the model is that the classical binding
energies predicted for nuclei are too big. One possible strategy to reduce the
binding energies is to look for modifications of the model that posses self-
dual sectors. There have been several attempts in the literature to confront
this problem either by adding higher-derivative terms adam-BPS-Skyrme , or by
adding extra fields sutcliffe2 ; sutcliffe3 .

In this work we will explore two extra possibilities. The first one is the
modification of the base space manifold. It is known that, the charge one
Skyrmion cannot be BPS if the space manifold is $\mathbb{R}^{3}$. This is
because the hypothetical BPS solution would be an isometry between the base
and target space manifolds and there are not isometries between
$\mathbb{R}^{3}$ and $\mathbb{S}^{3}$. If one considers the formulation of the
model in $\mathbb{S}^{3}$, the charge one Skyrmion is the identity map between
three-spheres, which is trivially an isometry and it is automatically BPS
Manton2 . The second possibility is to allow the coupling constants to be
“running couplings” or dielectric terms Adam1 . This allows again for charge
one BPS Skyrmions. In the present work we generalize these results to higher-
charge Skyrmions.

We will also study a gravitating Skyrme model. It is remarkable that, the
gravitating $B=1$ case it is actually solvable Canfora1 ; Canfora2  and self-
dual in the matter sector, provided that a cosmological constant is added to
the model. We will also extend this result to higher-charges.

This paper is organized as follows. In Sec. 2 we describe a Skyrme model with
BPS sectors in curved space. In Sec. 3, we analyze several dielectric Skyrme
models with arbitrary BPS sectors. In Sec. 4 we study the gravitating Skyrme
model with cosmological constant. Finally, Sec. 5 contains our conclusions and
further discussion.

##  2 The Skyrme model in curved spacetime

In this paper we will focus on static configurations, therefore we will mainly
be interested in three dimensional manifolds (spatial slices at fixed time of
the four dimensional spacetime manifold). The original Skyrme model consists
of two terms, a Dirichlet type term, quadratic in derivatives, and a Skyrme
term, quartic in derivatives. The static energy functional over an arbitrary
three-dimensional manifold $\mathcal{M}_{3}$, characterized by a metric
$g_{ij}$, has the following form

$E=\int
d^{3}x\sqrt{-g}\left(-\frac{f^{2}}{2}g^{ij}\text{tr}\,R_{i}R_{j}-\frac{1}{16e^{2}}g^{ij}g^{kl}\text{tr}[R_{i},R_{k}][R_{j},R_{l}]\right),$
(2.1)

where $R_{i}=\partial_{i}UU^{-1}$ are the $su(2)$-valued currents, $U\in
SU(2)$ contains the fundamental degrees of freedom and $e,f$ are constants. It
is very useful to define the strain tensor $D_{ij}$ as follows

$D_{ij}=-\frac{1}{2}\text{tr}R_{i}R_{j}.$ (2.2)

By using the eigenvalues of $D_{ij}$ the static energy functional can be cast
as

$E=\int
d^{3}x\sqrt{-g}\left(f\omega_{m}\pm\frac{1}{2e}|\epsilon_{mnl}|\omega_{n}\omega_{l}\right)^{2}\mp\frac{f}{e}6\int
d^{3}x\lambda_{1}\lambda_{2}\lambda_{3},$ (2.3)

where $\omega_{m}=e^{\,i}_{m}\lambda_{i}$, $e^{i}_{m}$ is the dreibein
verifying $g^{ij}=e^{i}_{m}e^{j}_{m}$. The energy functional (2.3) gives a
lower bound

$|E|\geq\frac{f}{e}6\int d^{3}x\lambda_{1}\lambda_{2}\lambda_{3},$ (2.4)

which is saturated when the first term in (2.3) vanishes. This gives the
following set of BPS equations

$f\omega_{m}\pm\frac{1}{2e}|\epsilon_{mnl}|\omega_{n}\omega_{l}=0.$ (2.5)

Of course, in general it is not guaranteed that equations (2.5) have
nontrivial solutions. For example, as we have said, in $\mathbb{R}^{3}$, and
for $e,f$ constants, (2.5) are only verified for constant fields.

###  2.1 BPS structure

We will analyze under which circumstances the equations (2.5) can be
satisfied. Since we are interested in looking for a single positive answer to
the problem, we will restrict the possible solutions to only those that
respect a particular form. As we will see later, this choice not only works as
a pro-example, but also allows us to find trivial analytical solutions on each
topological sector. Let us consider, then, the following unit vector

$\hat{\mathbb{n}}_{R}=\frac{1}{1+|R|^{2}}\left(2\,\text{Re}\,R,2\,\text{Im}\,R,1-|R|^{2}\right),$
(2.6)

where $R$ is a rational function of $z$. We define the $SU(2)$ Skyrme field in
the rational map ansatz (Manton1 , Manton3 ) as follows

$U(r,z)=e^{i\xi(r)\hat{\mathbb{n}}_{R}\cdot\bar{\mathbb{\sigma}}},$ (2.7)

where $\bar{\mathbb{\sigma}}=(\sigma_{1},\sigma_{2},\sigma_{3})$ are the Pauli
matrices. The eigenvalues of the strain tensor (2.2) for the ansatz (2.7) are
given by

$\lambda_{1}^{2}=\xi^{\prime
2}(r),\,\,\,\lambda_{2}^{2}=\lambda_{3}^{2}=\frac{\sin^{2}\xi(r)}{r^{2}}\frac{\left(1+|z|^{2}\right)^{2}}{\left(1+|R|^{2}\right)^{2}}\mathinner{\\!\left\lvert\frac{dR}{dz}\right\rvert}^{2}.$
(2.8)

Taking into account the form of the eigenvalues of the strain tensor given in
(2.8) and the general form of the (possible) BPS equations, we consider the
following base space metric

$g_{ij}=\frac{4r_{0}^{2}}{\left(1+r^{2}\right)^{2}}\text{diag}\left(1,r^{2}\left(\frac{1+|z|^{2}}{1+|Q|^{2}}\right)^{2}\mathinner{\\!\left\lvert\frac{dQ}{dz}\right\rvert}^{2},r^{2}\left(\frac{1+|z|^{2}}{1+|Q|^{2}}\right)^{2}\mathinner{\\!\left\lvert\frac{dQ}{dz}\right\rvert}^{2}\sin^{2}\theta\right).$
(2.9)

where $z=\tan\left(\frac{\theta}{2}\right)e^{i\phi}$ and $Q$ is a fixed
rational function of $z$. Note that for $Q(z)=z$, and after the change of
coordinates $r=\tan\psi/2$, the metric (2.9) takes the form of the round
metric in $\mathbb{S}^{3}$ of radius $r_{0}$

$ds^{2}=r_{0}^{2}\left(d\psi^{2}+\sin^{2}\psi(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\right).$ (2.10)

The BPS equations (2.5) with the metric (2.9) take the form

$\displaystyle
ef\xi^{\prime}(r)+\frac{\left(1+|Q|^{2}\right)^{2}}{\left(1+|R|^{2}\right)^{2}}\frac{|R^{\prime}|^{2}}{|Q^{\prime}|^{2}}\frac{1+r^{2}}{r_{0}r^{2}}\sin^{2}(\xi(r))$
$\displaystyle=$ $\displaystyle 0$ (2.11)
$\displaystyle\frac{1+|Q|^{2}}{1+|R|^{2}}\frac{|R^{\prime}|}{|Q^{\prime}|}\left(\frac{1}{r_{0}}\xi^{\prime}(r)+2\frac{ef}{1+r^{2}}\right)$
$\displaystyle=$ $\displaystyle 0.$ (2.12)

If $ef=1/r_{0}$ these equations have a common solution

$\xi(r)=-2\arctan r+\pi,\quad R(z)=Q(z),$ (2.13)

which interpolates between $\pi$ and $0$ for $r\in(0,\infty)$. In terms of the
coordinate $\psi\in[0,\pi]$ the solution is

$\xi(\psi)=\pi-\psi,$ (2.14)

which represents the identity map. This solution was previously found for
$B=1$ in $\mathbb{S}^{3}$ in Manton2  and in terms of the dielectric functions
in Adam1  (in the latter, this solution corresponds to the choice $\alpha=1/2$
of that paper). The topological charge of the solution is actually determined
by the choice of $R(z)$ vie the usual relation

$B=-\frac{1}{2\pi^{2}}\int\frac{\xi^{\prime}(r)\sin^{2}\xi(r)}{\left(1+|r|^{2}\right)^{2}}\mathinner{\\!\left\lvert\frac{dR}{dz}\right\rvert}^{2}2id\bar{z}dzdr=N,$
(2.15)

where $N$ is the degree of the map $R$. As a consequence, not only the $B=1$
sector can be made BPS by choosing the round $\mathbb{S}^{3}$ to be the base
space manifold, but any choice of (2.9) with $\text{deg}Q=N$ brings the $B=N$
sector to self-duality. Moreover, this BPS solution is given by the rational
map ansatz defined by the same rational function defining the metric. This is
behind the choice of the metric (2.9), which is the unique metric that allows
for exact BPS solutions in the form of the rational map ansatz. In addition,
it is straightforward to prove that (2.11) and (2.12) imply the Euler-Lagrange
equations. Another feature of the solutions (2.13) is that the length scale
determined by the Syrmion solutions, $L=1/ef$ must equal the radius $r_{0}$ of
the base space manifold is order for the solutions to be BPS. It is worth
noting that for any choice of metric on the form (2.9) the radial profile of
the solutions in the self-dual sector is of the form given in (2.13), i.e. it
is an identity map in the $\psi$ coordinate regardless of the topological
charge, while the angular part is given by $Q(z)$. In principle, there may be
other metric choices, not related to the rational map ansatz providing BPS
sectors as well. The strategy for finding them follows the lines described
above: fix another field ansatz, compute the eigenvalues for the strain tensor
and choose $g_{ij}$ such that the BPS equations as satisfied simultaneously.

##  3 Dielectric Skyrmions with higher-charge BPS sectors

We have considered so far, how modifications of the space manifold may
generate certain BPS sectors (of any baryonic charge). That is, a judicious
modification of the space metric in which the Skyrme model is formulated
allows for first order equations saturating a Bogomolny bound. We discuss now
if some “mild” modification of the Skyrme model may lead to similar results,
but now keeping the base space manifold to be $\mathbb{R}^{3}$. One
possibility is the addition of impurities (background space-dependent fields).
This was already considered in Naya  (see also tong ; BPS-imp-susy ; solvable-
imp ; vortex-imp-2 ; vortex-imp-3  for other BPS models with imourities), but
we will not pursue this line. In Adam1 , the authors consider a dielectric
Skyrme model, that consists of the original Skyrme model (quadratic plus
quartic terms), multiplied by some dielectric functions. These dielectric
functions depend on the field variables and therefore, they can be viewed as a
sort of modification of the target space metric, even though these dielectric
functions allow in general for more freedom than target-space metric
modifications. They found that, with such a modification, it is possible to
bring the $B=1$ Skyrmion to the BPS sector. We will show here that it is
indeed possible to generate BPS sectors of arbitrary charge (but only one at a
time) by modifying properly the dielectric functions. Let us consider the
following action

$E=\int
d^{3}x\left(-\frac{1}{2}G_{t}^{ij}(\pi^{a})\text{tr}\,R_{i}R_{j}-\frac{1}{16}G_{t}^{ij}(\pi^{a})G_{t}^{kl}(\pi^{a})\text{tr}[R_{i},R_{j}][R_{k},R_{l}]\right),$
(3.1)

where $G_{t}^{ij}(\pi^{a})$ are functions of the target space variables,
generically denoted by $\pi^{a}$, defined as follows

$G_{t}^{ij}(\pi^{a})=\text{diag}(f_{1}^{2}(\pi^{a}),f_{2}^{2}(\pi^{a}),f_{3}^{2}(\pi^{a})).$
(3.2)

Note that we have reabsorbed the constant couplings $e,f$ in the definitions
of the functions $f_{i}$. By completing the square, the energy (3.1) we can be
expanded as follows

$\displaystyle E$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{3}}\left(f_{1}\lambda_{1}\pm
f_{2}f_{3}\lambda_{2}\lambda_{3}\right)^{2}+\left(f_{2}\lambda_{2}\pm
f_{1}f_{3}\lambda_{1}\lambda_{3}\right)^{2}+\left(f_{3}\lambda_{3}\pm
f_{1}f_{2}\lambda_{1}\lambda_{2}\right)^{2}$ (3.3) $\displaystyle\mp
6\int_{\mathbb{R}^{3}}f_{1}f_{2}f_{3}\lambda_{1}\lambda_{2}\lambda_{3},$

therefore

$|E|\geq
6\int_{\mathbb{R}^{3}}f_{1}f_{2}f_{3}\lambda_{1}\lambda_{2}\lambda_{3}.$ (3.4)

The bound is saturated if and only if

$(f_{1}\lambda_{1})^{2}=(f_{2}\lambda_{2})^{2}=(f_{3}\lambda_{3})^{2}=1.$
(3.5)

It is well-known that, in the standard case this bound cannot be attained for
any nontrivial configuration and, as we have said, it is also known that there
is a choice of dielectric functions Adam1  that makes the $B=1$ sector BPS.
Let us explore the rest of sectors. For the sake of concreteness let us take
the rational map ansatz (2.7) whose target space coordinates $\pi^{a}$ are
expanded in terms of $\xi(r)$ and $R(z)$. The relations (3.5) become

$\displaystyle f_{1}^{2}(\xi,R)\xi^{\prime}(r)^{2}$ $\displaystyle=$
$\displaystyle 1$ (3.6) $\displaystyle
f_{2}^{2}(\xi,R)\frac{\sin^{2}\xi(r)}{r^{2}}\frac{\left(1+|z|^{2}\right)^{2}}{\left(1+|R|^{2}\right)^{2}}\mathinner{\\!\left\lvert\frac{dR}{dz}\right\rvert}^{2}$
$\displaystyle=$ $\displaystyle 1$ (3.7) $\displaystyle\
f_{3}^{2}(\xi,R)\frac{\sin^{2}\xi(r)}{r^{2}}\frac{\left(1+|z|^{2}\right)^{2}}{\left(1+|R|^{2}\right)^{2}}\mathinner{\\!\left\lvert\frac{dR}{dz}\right\rvert}^{2}$
$\displaystyle=$ $\displaystyle 1.\ $ (3.8)

Note that (3.6)-(3.8) are in fact very similar to (2.5), but there is a
crucial difference: in the first case the modification is in terms of
spacetime coordinates (dreibeins to be precise), while in the latter, the
functions $f_{i}$ depend on the target space coordinates, that is, in the
fields themselves. These relations imply the following constraints on the
functions $f_{i}$

$\displaystyle f_{1}(\xi,R)$ $\displaystyle=$ $\displaystyle f_{1}(\xi),$
(3.9) $\displaystyle f_{2}(\xi,R)$ $\displaystyle=$ $\displaystyle
f_{3}(\xi,R).$ (3.10)

In addition, if we assume that $f_{2}(\xi,R)$ factorizes, i.e.
$f_{2}(\xi,R)=f_{2}(\xi)F(R,\bar{R})$ we have

$\displaystyle f_{2}^{2}(\xi)\frac{\sin^{2}\xi(r)}{r^{2}}$ $\displaystyle=$
$\displaystyle C,$ (3.11) $\displaystyle
F(R,\bar{R})\frac{\left(1+|z|^{2}\right)^{2}}{\left(1+|R|^{2}\right)^{2}}\mathinner{\\!\left\lvert\frac{dR}{dz}\right\rvert}^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{C},$ (3.12)

where $C$ is a constant that we will choose to one for simplicity. Let us
analyze first the angular equation (3.12). A judicious choice $F$
automatically determines $R$ and define the topological charge of the
solution. For example

$\displaystyle F(R,\bar{R})=1$ $\displaystyle\Rightarrow$ $\displaystyle
R(z)=z,$ (3.13) $\displaystyle
F(R,\bar{R})=\frac{(1+|R|^{2})^{2}}{4|R|(1+|R|)^{2}}$
$\displaystyle\Rightarrow$ $\displaystyle R(z)=z^{2},$ (3.14)

which correspond respectively to a charge 1 and charge 2 BPS solutions.
Explicit expressions for $F(R,\bar{R})$ of topological charges 3 and 4 are not
specially illuminating, since the involve inversions of 3 and 4 degree
polynomials. For even higher topological charges an explicit expression is not
guaranteed. In general, one can compute the corresponding $F(R,\bar{R})$ of
any topological charge and with a prescribed symmetry in the following way.
First, we choose a rational map $R=S(z)$ of any degree and invert the
equation, $z=S^{-1}(R)$. Then, we choose $F(R,\bar{R})$ such that (3.12) is
satisfied. The obvious choice is

$F(R,\bar{R})=\frac{(1+R\bar{R})^{2}}{(1+S^{-1}(R)S^{-1}(\bar{R}))^{2}}|R^{\prime}(S^{-1}(R))|^{-2}.$
(3.15)

Note that in this way, $F(R,\bar{R})$ is entirely written in terms of target
space variables and satisfies (3.12). It still remains to determine the
functions $f_{1}(\xi)$ and $f_{2}(\xi)$. In principle we can choose almost
freely one of the functions, let us say $f_{1}(\xi)$. Then, by consistency,
$f_{2}(\xi)$ must be chosen in such a way that (3.11) is verified by the
solution of (3.6). The consistency condition (remember that we choose $C=1$)
between both functions may be written symbolically as follows

$\int f_{1}(\xi)d\xi=f_{2}(\xi)\sin\xi.$ (3.16)

For example, we can reproduce the solution of Sec. 2 with the choice

$f_{1}^{2}(\xi)=\frac{1}{(\cos\xi-1)^{2}}\Rightarrow f_{2}(\xi)=f_{1}(\xi)$
(3.17)

Actually, once a choice of the profile functions has been made, it generates a
solution in any BPS sector provided that we choose properly the angular
function $F(R,\bar{R})$. It other words, it is the choice of the function
$F(R,\bar{R})$ that determines which topological sector is BPS. It is clear
that all the solutions constructed in this way saturate the energy bound (3.5)
and it is straightforward (but tedious) to prove that the also verify the
Euler-Lagrange equations from (3.1). It must be emphasized that the choice
(3.17) is far from unique. If $f_{1}(\xi)$ is chosen in such a way that (3.6)
satisfies the correct boundary conditions, then $f_{2}(\xi)$ can be computed
formally from (3.16).

##  4 The gravitating Skyrme model

As we have seen, for a fixed spacetime manifold the Skyrme model admits one
BPS sector of arbitrary degree. Let us consider now the coupling to gravity.
In this case we have

$\mathcal{L}=\frac{1}{2\kappa}\int\sqrt{-g}\left(\mathcal{R}-2\Lambda\right)+\int\sqrt{-g}\left(\frac{f^{2}}{2}\mathcal{L}_{2}+\frac{1}{16e^{2}}\mathcal{L}_{4}\right),$
(4.1)

where $\Lambda$ is a cosmological constant. We will consider the following
metric ansatz

$ds^{2}=-dt^{2}+g_{ij}dx^{i}dx^{j},$ (4.2)

where $g_{ij}$ is the base space metric and $x^{i}=(r,\theta,\phi)$. The
solutions of the form (2.7) with $\xi(r)$ defined by (2.13) still saturate the
energy bound (2.4) for a metric of the form (2.9). But now, due to the
coupling to gravity, the metric is not fixed, but it must satisfy the Einstein
equations

$G_{\mu\nu}+\Lambda g_{\mu\nu}=\kappa T_{\mu\nu},$ (4.3)

where $G_{\mu\nu}$ is the Einstein tensor and $T_{\mu\nu}$ is the energy-
momentum tensor. Let us consider then, the family of 3-manifolds,
$\mathcal{M}_{3}(Q)$ with metric (2.9) (labeled by the rational function $Q$).
The volume and the curvature scalar can be computed easily

$\text{Vol}(\mathcal{M}_{3}(Q))=2\pi
r_{0}^{2}\text{deg}Q,\quad\mathcal{R}(\mathcal{M}_{3}(Q))=\frac{6}{r_{0}^{2}}.$
(4.4)

Now, if we assume that a metric of the form (2.9) is a solution of the
Einstein equations then, the Einstein tensor has a simple form

$\displaystyle G_{0}^{\,0}$ $\displaystyle=$
$\displaystyle-\frac{3}{r_{0}^{2}},$ (4.5) $\displaystyle G_{i}^{\,i}$
$\displaystyle=$ $\displaystyle-\frac{1}{r_{0}^{2}}.$ (4.6)

From Sec. 2, we know that the BPS solutions in the metric (2.9) satisfy (2.11)
and (2.12), provided that they are given by the rational map ansatz. The
energy-momentum tensor in the rational map ansatz can be expanded as follows

$\displaystyle T_{t}^{\,t}$ $\displaystyle=$
$\displaystyle\Delta_{1}-\xi^{\prime 2}(r)\Delta_{2},$ (4.7) $\displaystyle
T_{r}^{\,r}$ $\displaystyle=$ $\displaystyle\Delta_{1}+\xi^{\prime
2}(r)\Delta_{2},$ (4.8) $\displaystyle T_{\theta}^{\,\theta}$ $\displaystyle=$
$\displaystyle\frac{(1+r^{2})^{2}}{16r_{0}^{4}}\left(\frac{(1+r^{2})^{2}\sin^{4}\xi(r)}{e^{2}r^{4}}-4f^{2}r_{0}^{2}\xi^{\prime
2}(r)\right),$ (4.9) $\displaystyle T_{\phi}^{\,\phi}$ $\displaystyle=$
$\displaystyle T_{\theta}^{\,\theta},$ (4.10)

where

$\displaystyle\Delta_{1}$ $\displaystyle=$
$\displaystyle\frac{(1+r^{2})^{2}}{32e^{4}r^{4}r_{0}^{4}}\sin^{2}\xi(r)\left((1+r^{2})^{2}\cos
2\xi(r)-1-r^{2}(2+r^{2}+16e^{2}f^{2}f_{0}^{2})\right),$ (4.11)
$\displaystyle\Delta_{2}$ $\displaystyle=$
$\displaystyle\frac{(1+r^{2})^{2}}{8e^{2}r^{2}r_{0}^{4}})\left((1+r^{2})^{2}\sin^{2}\xi(r)+2e^{2}f^{2}r^{2}r_{0}^{2}\right).$
(4.12)

It should be noted that the energy-momentum tensor does not depend on the
rational function $R$, but only on the profile function. Our hypothesis is
that the BPS solutions already found for a fixed background metric (2.13) are
still solutions in the gravitating case. Therefore, one should restrict the
energy-momentum tensor to the BPS sector. The on-shell energy-momentum tensor
has the form

$\displaystyle T_{0}^{\,0}$ $\displaystyle=$
$\displaystyle-3\frac{1+f^{2}e^{2}r_{0}^{4}}{e^{2}r_{0}^{4}},$ (4.13)
$\displaystyle T_{i}^{\,i}$ $\displaystyle=$
$\displaystyle-\frac{1-e^{2}f^{2}r_{0}^{2}}{e^{2}r_{0}^{4}}.$ (4.14)

Finally, taking into account (4.5)- (4.6) and (4.13)-(4.14), Einstein
equations (4.3) are identically satisfied if the following constraints for the
radius of the base space and the cosmological constant hold

$r_{0}^{2}=\frac{2\kappa}{e^{2}(1-\kappa
f^{2})},\quad\Lambda=\frac{3}{4\kappa}e^{2}\left(\kappa f^{2}-1\right)^{2}.$
(4.15)

This is the consistency condition we were looking for. As long as (4.15)
holds, field equations and Einstein equations are satisfied simultaneously for
a rational function $R$ of any degree. Similar constraints were obtained for
the $B=1$ case in Canfora1 ; Canfora2 . Moreover, since the condition
$ef=1/r_{0}$ must be met for the BPS equations to be satisfied we have

$\Lambda=e^{2}f^{2},\quad\kappa=\frac{1}{3f^{2}}.$ (4.16)

On constant time hypersurfaces, the topological charge of a BPS soliton is
well-defined and determined by the degree of $R$ provided that (4.16) are
satisfied. Therefore, the model (4.1) has solutions of arbitrary degree
consistent with Einstein equations. It should be noted that for any degree,
the radial profile is given by (2.13) with an angular dependence encoded in
the rational map ansatz. It is worth noting that the cosmological constant
cannot vanish due to the constraints (4.16). A vanishing $\Lambda$ implies
either a vanishing Dirichlet term or a divergent Skyrme term. Moreover, the
base space manifold must be compact ($r_{0}$ finite) since the condition
$ef=1/r_{0}$ must hold in order to have BPS solutions. In addition, the size
of the Skyrmion $L=r_{0}=1/\sqrt{\Lambda}$ has to cover the whole base space
in order the solution to be BPS. This feature was also observed in Manton2
for the $B=1$ Skyrmion in $\mathbb{S}^{3}$ and in Gudnason  for the sextic
(BPS) Skyrme model.

##  5 Conclusions

We have discussed two modifications of the Skyrme model allowing for self-dual
sectors. First, we have proposed a family of curved space manifolds labeled by
rational maps. A particular choice of such a map automatically fixes one (and
only one) BPS sector of topological charge $N$, being $N$ the degree of the
rational map. Moreover, for all possible choices of rational maps, the radial
part of the solution takes the same form, and corresponds to the identity map
after a suitable change of coordinates. The angular part of the solutions
corresponds exactly to the rational map chosen to build the space metric in
the first place.

In the second part of the paper we have discussed a different deformation of
the Skyrme model. This time, instead of curving the space manifold we have
added a set of dielectric functions. This could be interpreted as a sort of
modification of the target space metric. We have demonstrated that, a
particular choice of these field dependent functions brings one topological
sector to self-duality. We have also provided an explicit example, based on
the rational map ansatz, where the BPS sector in exactly solvable.

Regarding possible applications to nuclear physics, one should look for near-
BPS models rather than strict BPS models. In our case, it is not very
difficult to convince oneself that small deformations of the dielectric
functions will move the corresponding sector from BPS to near-BPS. A related
question concerns the binding energies of those sectors that were not BPS in
the original model. The huge freedom in the choice of $F$ and $f_{i}$ makes
very difficult to answer this question in general. At least, for one
particular choice Adam1  the binding energies can be reduced to a $0.2\%$.
Whether or not it is possible to obtain realistic binding energies for other
choices of the dielectric functions for a range of baryonic charges will
depend of course on the details of the dielectric functions, by we can say a
priori that, within this family of models, at least one BPS sector (of any
charge) will have arbitrarily small binding energy, depending on how the near-
BPS deformation is done. On the other hand, a physical motivation for the
dielectric functions, whether they are obtained from higher-order terms in the
effective field theory perspective or as in-medium coupling constants, may
reduce this freedom. We think that these issues deserve further investigation.

It should be noted that the models provided in Secs. 2 and 3 have only one
self-dual sector, i.e. they do not correspond to a standard BPS model (where
infinite topological sectors give rise to infinite self-dual sectors).
However, the gravitating case discussed in Sec. 4 seems to have this property.
The Einstein equations impose some constraints between the size of the
spacetime manifold and the cosmological constant, or more precisely a relation
between both of them and the couplings of the model. In this situation
Einstein equations are satisfied for a metric of the form (2.9) and the
rational map ansatz automatically satisfies the BPS equations in all sectors.

Another interesting feature that deserves further study is the possible
existence of a moduli space. As we have mentioned, in the non gravitating
case, a choice of rational map $R(z)$ for the space metric, determines the BPS
sector and the exact solution in the BPS sector is defined by the rational map
itself. This means that, once $R(z)$ has been chosen, say of degree $N$, any
other rational map of the same (or different degree) is not a BPS solution.
This implies in particular that there is no moduli space associated to that
BPS sector. The situation in the gravitating case is different. In Sec. 4 we
have shown that, once a cosmological constant has been added, the gravitating
model has infinitely many BPS sectors (of any topological charge) whose BPS
solutions are given by rational maps and metric solutions given by (2.9).
Since the BPS sector is determined by the degree of the rational map, a
trajectory on the moduli space corresponds to the variation of the
coefficients of the rational map (preserving the degree), which in turn,
implies a variation of the metric of the base space. A detailed analysis of
these moduli spaces is left for a future work.

Finally, it would be also interesting to explore the possible supersymmetric
structure of the gravitating model presented here. In three spacetime
dimensions there are BPS Skyrme models, that is, Skyrme models with infinite
self-dual sectors. It is well-known that there is a relation between self-
duality and SUSY, and in fact, SUSY versions of these low dimensional models
have been found Queiruga1 ; Queiruga2 ; Queiruga3 ; Bolognesi  some years ago.
So far, the attempts of supersymetrization of the Skyrme model in four
dimensions were only partial Nepomechie ; Queiruga4 ; Gudnason1 . Actually,
the supersymmetric form of the Skyrme term alone, both in two and three
dimensions, requires the introduction of a constant (Queiruga2 ; Gudnason1 )
that depends on the couplings as (4.16), but it does not allow for a Dirichlet
term (although in two dimensions this constant can be promoted to a field
dependent potential preserving SUSY). The gravitating Skyrme model described
in Sec. 4 suggests that, there should be a SUGRA formulation of this model.
These problems are also left for a future investigation.

Acknowledgements.\- This work is supported by the Spanish Ministry
MCIU/AEI/FEDER grant (PGC2018-094626-B-C21), the Basque Government grant
(IT-979-16). The author is grateful to C. Adam and A. Wereszczynski for useful
comments.

## References

  * (1) T.H.R. Skyrme, “A nonlinear field theory,” Proc. R. Soc. Lond. 260 (1961) 127.
  * (2) T.H.R. Skyrme, “A unifield field theory of mesons and baryons,” Nucl. Phys. 31 (1962) 556.
  * (3) E. Witten, “Baryons in the 1/N expansion,” Nucl. Phys. B 160 (1979) 57.
  * (4) E. Witten, “Current algebra, baryons, and quark confinement,” Nucl. Phys. B 223 (1983) 433.
  * (5) G.S. Adkins, C.R. Nappi and Edward Witten, “Static Properties of Nucleons in the Skyrme Model,” Nucl.Phys.B 228 (1983) 552.
  * (6) C. Adam, J. Sanchez-Guillen and A. Wereszczynski, “A Skyrme-type proposal for baryonic matter,” Phys.Lett.B 691 (2010) 105-110, arXiv:1001.4544 [hep-th].
  * (7) P.M. Sutcliffe, “Skyrmions, instantons and holography,” JHEP 08 (2010) 019, arXiv:1003.0023 [hep-th].
  * (8) P.M. Sutcliffe, “Skyrmions in a truncated BPS theory,” JHEP 04 (2011) 045, arXiv:1101.2402 [hep-th].
  * (9) N.S. Manton and P.J. Ruback, “Skyrmions in Flat Space and Curved Space,” Phys.Lett.B 181 (1986) 137-140.
  * (10) S.B. Gudnason, M. Nitta and N. Sawado, “Gravitating BPS Skyrmions”, JHEP 12 (2015) 013, arXiv: 1510.08735.
  * (11) C. Adam, K. Oles and A. Wereszczynski, “The Dielectric Skyrme model,” Phys.Lett.B 807 (2020) 135560, arXiv: 2005.00018 [hep-th].
  * (12) N.S. Manton and P. Sutcliffe, “Topological Solitons,” (Cambridge University Press, Cambridge, 2007).
  * (13) C. J. Houghton, N.S. Manton and P.M. Sutcliffe, “Rational maps, monopoles and Skyrmions,” Nucl.Phys.B 510 (1998) 507-537, arXiv: hep-th/9705151 [hep-th].
  * (14) E. Ayon-Beato, F. Canfora and J. Zanelli, “Analytic self-gravitating Skyrmions, cosmological bounces and AdS wormholes,” Phys.Lett.B 752 (2016) 201-205 , arXiv: 1509.02659 [gr-qc].
  * (15) F. Canfora, A. Paliathanasis, T. Taves and J. Zanelli, “Cosmological Einstein-Skyrme solutions with nonvanishing topological charge,” Phys.Rev.D 95 (2017) 6, 065032, arXiv: 1703.04860 [hep-th].
  * (16) C. Naya and K. Oles, “Background fields and self-dual Skyrmions,” Phys.Rev.D 102 (2020) 2, 025007, arXiv: 2004.07069 [hep-th]
  * (17) D. Tong, K. Wong, Vortices and Impurities, JHEP 1401 (2014) 090.
  * (18) C. Adam, J. Queiruga and A. Wereszczynski, BPS soliton-impurity models and supersymmetry, JHEP 1907 (2019) 164.
  * (19) C. Adam, K. Oles, J. Queiruga, T. Romanczukiewicz and A. Wereszczynski, Solvable self-dual impurity models, JHEP 1907 (2019) 150.
  * (20) A. Cockburn, S. Krusch, A. Muhamed, Dynamics of vortices with magnetic impurities, J.Math.Phys. 58 (2017) 063509.
  * (21) J. Ashcroft, S. Krusch, Vortices and magnetic impurities, Phys.Rev.D 101 (2020) 2, 025004 ,arXiv:1808.07441 [hep-th].
  * (22) C. Adam, J.M. Queiruga, J. Sanchez-Guillen and A. Wereszczynski, “$N=1$ supersymmetric extension of the baby Skyrme model,” Phys.Rev. D84 (2011) 025008, arXiv:1105.1168 [hep-th].
  * (23) C. Adam, J.M. Queiruga, J. Sanchez-Guillen and A. Wereszczynski, “Extended supersymmetry and BPS solutions in baby Skyrme models,” JHEP 1305 (2013) 108, arXiv:1304.0774 [hep-th].
  * (24) J. M. Queiruga, “Baby Skyrme model and fermionic zero modes,” Phys.Rev. D94 (2016) no.6, 065022, arXiv:1606.02869 [hep-th].
  * (25) S. Bolognesi and W. Zakrzewski, “Baby Skyrme model, near-BPS approximations and supersymmetric extensions,” Phys.Rev. D91 (2015) no.4, 045034, arXiv:1407.3140 [hep-th].
  * (26) E. A. Bergshoeff, R.I. Nepomechie and H.J. Schnitzer, ”Supersymmetric Skyrmions in four-dimensions,” Nucl.Phys.B 249 (1985) 93-130.
  * (27) J.M. Queiruga, “Skyrme-like models and supersymmetry in 3+1 dimensions,” Phys.Rev.D 92 (2015) 10, 105012, arXiv:1508.06692 [hep-th].
  * (28) S.B. Gudnason, M. Nitta and S. Sasaki, “A supersymmetric Skyrme model,” JHEP 02 (2016) 074, arXiv:1512.07557 [hep-th].