# The cos$\phi_{R}$ asymmetry of dihadron production in double longitudinally
polarized SIDIS

Wei Yang  School of Physics, Southeast University, Nanjing 211189, China  Hui
Li  School of Physics, Southeast University, Nanjing 211189, China  Zhun Lu
<EMAIL_ADDRESS>School of Physics, Southeast University, Nanjing 211189,
China

###### Abstract

We present a study on the double longitudinal-spin asymmetry of dihadron
production in semi-inclusive deep inelastic scattering, in which the total
transverse momentum of the final-state hadron pairs is integrated out. In
particular, we investigate the origin of the $\cos\phi_{R}$ azimuthal
asymmetry for which we take into account the coupling of the helicity
distribution $g_{1}$ and the twist-3 dihadron fragmentation function
$\widetilde{D}^{\sphericalangle}$. We calculate the $s$-wave and $p$-wave
interference term $\widetilde{D}^{\sphericalangle}_{ot}$ in a spectator model.
We estimate the $\cos\phi_{R}$ asymmetry at the kinematics of COMPASS which is
collecting data on dihadron production in polarized deep inelastic scattering.
The prediction of the same asymmetry at JLab 12GeV and a future EIC are also
presented. Our study indicates that measuring the $\cos\phi_{R}$ asymmetry in
SIDIS may be a ideal way to probe the dihadron fragmentation function
$\widetilde{D}^{\sphericalangle}$.

##  I Introduction

The spin and azimuthal asymmetries in semi-inclusive deep inelastic scattering
(SIDIS) process Mulders:1995dh ; Kotzinian:1994dv ; Boer:1997nt ; Goeke:2005hb
; Bacchetta:2006tn  have been recognized as useful tools for exploring both
the partonic structure of the nucleon and the hadronization mechanism of
hadrons, which are among the main tasks in QCD and hadronic physics. According
to factorization Ji:2004wu , the cross-section of SIDIS process can be
expressed as the convolution of the parton distribution function, the
fragmentation function, and the hard scattering factor. Distribution and
fragmentation functions are important nonperturbative quantities encoding the
internal parton and spin structure of nucleon as well as the fragmentation
mechanism.

Recently, a lot of attention has also been paid on the higher-twist
contributions Efremov:2002ut ; Gamberg:2003pz ; Avakian:2007mv ; Mao:2012dk ;
Mao:2014aoa ; Mao:2014fma ; Song:2010pf ; Wei:2016far ; Yang:2016mxl ;
Lorce:2016ugb ; Bastami:2018xqd  in SIDIS. Although the rigorous proof on
factorization at the twist-3 level in SIDIS has not been achieved
Bacchetta:2008xw ; Gamberg:2006ru , it is shown Bacchetta:2006tn ;
Bacchetta:2003vn  that those effects are related to the twist-3 distributions
and fragmentation functions based on the tree-level calculation: the spin or
azimuthal dependent structure function can be expressed as the convolution of
the twist-3 distribution/fragmentation functions and the twist-2
fragmenation/distribution functions. Moreover, there are existing and ongoing
experimental measurements on the azimuthal asymmetries of single-hadron and
hadron pair (dihadron from a single jet) production at higher-twist level in
polarized SIDIS by COMPASS, HERMES, and CLAS at JLab, not only for single
hadron production Airapetian:2006rx ; Airapetian:2005jc ; Alekseev:2010dm ;
Aghasyan:2011ha ; Gohn:2014zbz ; Airapetian:2019mov , but also for dihadron
production S.Sirtl . For the later case, the azimuthal asymmetry at twist-2
level has been measured and was applied to extract Bacchetta:2011ip ;
Bacchetta:2012ty ; Radici:2015mwa ; Radici:2018iag  transversity from SIDIS
data Airapetian:2008sk ; Adolph:2012nw ; Adolph:2014fjw  based on the coupling
of transversity and the chiral-odd dihadron fragmentation function
$H_{1}^{\sphericalangle}$ Collins:1993kq ; Jaffe:1997hf ; Bacchetta:2002ux ;
Boglione:1999pz ; Boer:1999mm ; Bianconi:1999cd ; Radici:2001na ;
Courtoy:2012ry  in the collinear factorization.

In this work, we study the azimuthal asymmetry of dihadron production in
double longitudinally polarized SIDIS
$l^{\rightarrow}+p^{\rightarrow}\rightarrow l+h_{1}+h_{2}+X$ in the case the
total transverse momentum of the dihadron is integrated out. As shown in Ref.
Bacchetta:2003vn , in this particular process, there are two twist-3 terms
that might give rise to the asymmetry with a $\cos\phi_{R}$ modulation. The
first one is the coupling of the T-odd twist-3 distribution $e_{L}(x)$ and the
twist-2 dihadron $H_{1}^{\sphericalangle}$, while the second one is the
combination of the helicity distribution $g_{1}(x)$ and the twist-3 dihadron
fragmentation function (DiFF) $\tilde{D}^{\sphericalangle}$ originating from
quark-gluon-quark correlation. Here the symbol $\sphericalangle$ denotes that
the corresponding DiFF is the interference fragmentation function. However, if
the time reversal invariance is imposed and the gauge link is the only source
of a T-odd distribution, $e_{L}(x)$ should vanish Goeke:2005hb  and the
$e_{L}(x)\,H_{1}^{\sphericalangle}$ term will not contribute to the
$\cos\phi_{R}$ asymmetry. Based on this consideration, only one term, the
$g_{1}(x)\,\tilde{D}^{\sphericalangle}$ term, should give rise to the
asymmetry. This is different from the single-hadron production in SIDIS in
which usually several twist-3 terms contribute to one observable. Thus,
investigating the $\cos\phi_{R}$ asymmetry in double polarized SIDIS provides
an opportunity to study the unknown twist-3 DiFF $\tilde{D}^{\sphericalangle}$
in a less ambiguous way.

The remained content of the paper is organized as follows. In Sec. II, we will
briefly review the theoretical framework of the $\cos\phi_{R}$ asymmetry of
dihadron production in doubly polarized SIDIS. In Sec. III, we use a spectator
model to calculate the twist-3 dihadron fragmentation function
$\widetilde{D}^{\sphericalangle}$. As there is no measurement on
$\widetilde{D}^{\sphericalangle}$, model calculation is important for
accessing information on this unknown DiFF. In Sec. IV, we numerical estimate
the DIFF $\widetilde{D}^{\sphericalangle}$ and the $A^{\cos\phi_{R}}_{LL}$
double longitudinal-spin asymmetry at the kinematics of COMPASS, JLab 12 GeV
as well as EIC. Finally, in Sec. V, we provide conclusion for this paper.

##  II Formalism of $\cos\phi_{R}$ asymmetry in dihadron production in SIDIS

We consider the following dihadron SIDIS process,

$\displaystyle l^{\rightarrow}(\ell)+N^{\rightarrow}(P)\longrightarrow
l(\ell^{\prime})+h_{1}(P_{1})+h_{2}(P_{2})+X,$ (1)

in which a longitudinally polarized lepton collides on the longitudinally
polarized target proton $N$ via the exchange of a virtual photon. Here the
arrow $\rightarrow$ denotes the longitudinally polarization of the beam or the
proton target. The corresponding 4-momenta are given in parenthesis in the
above formula, then the virtual photon has the momentum
$q=\ell-\ell^{\prime}$. $P$ is the momentum of the target with mass $M$. In
this process, the final-state quark with momentum $k=p+q$ then fragments into
two final-state hadrons, $h^{+}$ and $h^{-}$, plus unobserved state X. The
momenta of the pair are denoted by $P_{1}$, $P_{2}$, respectively. We present
the following kinematical variables that are necessary to describe the
differential cross section and express the DiFFs,

$\displaystyle x=\frac{k^{+}}{P^{+}},\qquad y=\frac{P\cdot q}{P\cdot l},\qquad
z_{i}=\frac{P_{i}^{-}}{k^{-}},$ (2) $\displaystyle
z=\frac{P_{h}^{-}}{k^{-}}=z_{1}+z_{2},\qquad Q^{2}=-q^{2},\qquad s=(P+l)^{2},$
(3) $\displaystyle P_{h}=P_{1}+P_{2},\qquad R=(P_{1}-P_{2})/2,\qquad
M_{h}=\sqrt{P_{h}^{2}}.$ (4)

Furthermore, the 4-vectors are given in terms of the light-cone coordinates
$a^{\mu}=(a^{+},a^{-},a_{T})$, where the light-cone components are defined as
$a^{\pm}=(a^{0}\pm a^{3})/\sqrt{2}$, and $a_{T}$ is a bidimensional vector of
the transverse component. Therefore, $x$ represents the lightcone momentum
fraction of the initial quark, $z_{i}$ is the lightcone momentum fraction of
hadron $h_{i}$ found in the fragmented quark. Finally, $M_{h}$, $P_{h}$ and
$R$ are the invariant mass, the total momentum and the relative momentum of
the hadron pair.

We work in the lab frame, in which the momenta $P_{h}^{\mu}$, $k^{\mu}$ and
$R^{\mu}$ can be decomposed to Bacchetta:2006un ; Boer:2003cm

$\displaystyle P_{h}^{\mu}$ $\displaystyle=$
$\displaystyle\left[P^{-}_{h},\frac{M_{h}^{2}}{2P^{-}_{h}},\vec{0}\right],$
(5) $\displaystyle k^{\mu}$ $\displaystyle=$
$\displaystyle\left[\frac{P_{h}^{-}}{z},\frac{z(k^{2}+\vec{k}_{T}^{2})}{2P_{h}^{-}},\vec{k}_{T}\right],$
(6) $\displaystyle R^{\mu}$ $\displaystyle=$
$\displaystyle\left[\frac{|\vec{R}|P^{-}_{h}}{M_{h}}\cos\theta,-\frac{|\vec{R}|M_{h}}{2P^{-}_{h}}\cos\theta,|\vec{R}|\sin\theta\cos\phi_{R},|\vec{R}|\sin\theta\sin\phi_{R}\right].$
(7) $\displaystyle=$
$\displaystyle\left[\frac{|\vec{R}|P^{-}_{h}}{M_{h}}\cos\theta,-\frac{|\vec{R}|M_{h}}{2P^{-}_{h}}\cos\theta,\vec{R}_{T}^{x},\vec{R}_{T}^{y}\right].$
(8)

Here, as shown in Fig. 1, $\phi_{R}$ is the angle between the lepton plane and
the dihadron plane defined as

$\displaystyle\cos{\phi_{R}}=\frac{\hat{q}\times\vec{l}}{|\hat{q}\times\vec{l}|}\cdot\frac{\hat{q}\times\vec{R_{T}}}{|\hat{q}\times\vec{R_{T}}|},\qquad\sin{\phi_{R}}=\frac{\vec{l}\times\vec{R_{T}}\cdot\hat{q}}{|\hat{q}\times\vec{l}||\hat{q}\times\vec{R_{T}}|}.$
(9)

where $\hat{q}=\vec{q}/|\vec{q}|$ and $R_{T}$ is the component of $R$
perpendicular to $P_{h}$. The angle $\theta$ is polar angle between the
direction of $P_{1}$ in the center of mass frame of the hadron pair and the
direction of $P_{h}$ in the lab frame.

There are several useful expressions for the scalar products as follows

$\displaystyle P_{h}\cdot R=0,$ (10) $\displaystyle P_{h}\cdot
k=\frac{M_{h}^{2}}{2z}+z\frac{k^{2}+|\vec{k}_{T}|^{2}}{2},$ (11)
$\displaystyle R\cdot
k=\left(\frac{M_{h}}{2z}-z\frac{k^{2}+|\vec{k}_{T}|^{2}}{2M_{h}}\right)|\vec{R}|\cos\theta-\vec{k}_{T}\cdot\vec{R}_{T}.$
(12)

In particular, there is a relation between the $\vec{R}$ and $M_{h}$.

$\displaystyle|\vec{R}|=\sqrt{\frac{M^{2}_{h}}{4}-m_{h}^{2}}.$ (13)

Figure 1: Sketch of the dihadron production in SIDIS process in the lab
framel, including the relevant azimuthal angles. The nucleon is assumed to be
longitudinally polarized either along or against the direction of the incoming
lepton.

We will focus on the case the lepton beam and proton target are both
longitudinally polarized as well as the case the total transverse momentum are
integrated out. For convenience, in the following formula we introduce the
subscript $U$, $L$ for the cross section $\sigma_{XY}$ to indicate the
unpolarized or longitudinally polarized states. Besides, the polarization with
respect to the beam and the target are indicated by the first label $X$ and
the second label $Y$ of $\sigma_{XY}$, respectively. We can then summarize the
unpolarized and double polarized cross sections as follows Ref.
Bacchetta:2003vn :

$\displaystyle\frac{d^{6}\\!\sigma_{UU}}{d\cos\theta\;dM_{h}^{2}\;d\phi_{R}\;dz\;dx\;dy}$
$\displaystyle=\frac{\alpha^{2}}{Q^{2}y}\,\left(1-y+\frac{y^{2}}{2}\right)\sum_{a}e_{a}^{2}f_{1}^{a}(x)\,D_{1}^{a}\bigl{(}z,M_{h}^{2},\cos\theta\bigr{)},$
(14)
$\displaystyle\frac{d^{6}\\!\sigma_{LL}}{d\cos\theta\;dM_{h}^{2}\;d\phi_{R}\;dz\;dx\;dy}$
$\displaystyle=\frac{\alpha^{2}}{Q^{2}y}\,S_{L}\sum_{a}e_{a}^{2}\,y(\frac{y}{2}-1)g^{a}_{1}(x)D^{a}_{1}(z,M_{h}^{2},\cos\theta)$
$\displaystyle+\frac{\alpha^{2}}{Q^{2}y}\,S_{L}\sum_{a}e_{a}^{2}\,2y\sqrt{2-y}{M\over
Q}\frac{|\bm{R}|}{M_{h}}\,\cos\phi_{R}$
$\displaystyle\times\left[xe^{a}_{L}(x)H_{1}^{\sphericalangle,a}\bigl{(}z,M_{h}^{2},\cos\theta\bigr{)}-\frac{M_{h}}{Mz}g^{a}_{1}(x)\widetilde{D}^{\sphericalangle,a}\bigl{(}z,M_{h}^{2},\cos\theta\bigr{)}\right].$
(15)

In Eq. (14), $f_{1}^{a}(x)$ and
$D_{1}^{a}\bigl{(}z,M_{h}^{2},\cos\theta\bigr{)}$ denote the unpolarized PDF
and unpolarized DiFF for flavor $a$. The first line in Eq.( II) represents the
leading twist contribution, while the second line denotes the twist-3 terms
which contributes to the $\cos\phi_{R}$ asymmetry. In details, there are two
individual contributions that might give rise to the asymmetry. The first one
is $e_{L}(x)H_{1}^{\sphericalangle}$, in which $e_{L}(x)=\int
d^{2}k_{T}e_{L}(x,k_{T}^{2})$ is a T-odd twist-3 distribution, and
$H_{1}^{\sphericalangle,a}$ is the twist-2 DiFF. The second one is
$g_{1}(x)\widetilde{D}^{\sphericalangle}$, where $g_{1}(x)$ is the twist-2
helicity distribution and $\widetilde{D}^{\sphericalangle}$ is the twist-3
DiFF, with the symbol $\sphericalangle$ denoting the interference nature of
the DiFF. However, as shown in Ref. Goeke:2005hb , if $e_{L}(x,k_{T}^{2})$
only receives contribution from the gauge-link, time-reversal invariance of
QCD implies the constraint $\int d^{2}k_{T}e_{L}(x,k_{T}^{2})=0$, therefore,
the contribution $e_{L}(x)H_{1}^{\sphericalangle}$ vanishes in the collinear
limit. Based on this observation, in this work we only need to consider the
$g_{1}(x)\widetilde{D}^{\sphericalangle}$ term.

The partial-wave analysis of the DiFFs $D_{1}$ and
$\widetilde{D}^{\sphericalangle}$ up to $p$-wave level yields Bacchetta:2003vn
:

$\displaystyle
D^{a}_{1}(z,\cos\theta,M_{h}^{2})=D^{a}_{1,oo}(z,M_{h}^{2})+D^{a}_{1,ol}(z,M_{h}^{2})\cos\theta+D^{a}_{1,ll}(z,M_{h}^{2})(3\cos^{2}\theta-1),$
(16) $\displaystyle H_{1}^{\sphericalangle
a}(z,\cos\theta,M_{h}^{2})=H_{1,ot}^{\sphericalangle
a}(z,M_{h}^{2})+H_{1,lt}^{\sphericalangle a}(z,M_{h}^{2})\cos\theta.$ (17)
$\displaystyle\widetilde{D}^{\sphericalangle}(z,\cos\theta,M_{h}^{2})=\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})+\widetilde{D}^{\sphericalangle}_{lt}(z,M_{h}^{2})\cos\theta.$
(18)

Here, $D^{a}_{1,oo}(z,M_{h}^{2})$ comes from the pure $s$-wave and $p$-wave
contributions, $D^{a}_{1,ol}$ and
$\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})$ arise from the
interference between a pair in $s$-wave and a pair in $p$-wave.

Following the similar arguments in Ref. Bacchetta:2006un , in this paper we
will not consider the $\cos\theta$-dependent terms in the expansion of DiFFs.
When integrating out the angular $\theta$ in the interval $[0,\pi]$ which is
our case, the $\cos\theta$-dependent terms should vanish. Therefore, we focus
on the functions $D^{a}_{1,oo}$ and $\widetilde{D}^{\sphericalangle}_{ot}$. In
this scenario, the $\cos\phi_{R}$ asymmetry of dihadron production
contribution to the double longitudinally polarized can be expressed as

$\displaystyle
A^{\cos\phi_{R}}_{LL}(x,z,M_{h}^{2})=-\frac{\sum_{a}e^{2}_{a}\frac{|\vec{R}|}{Q}\left[\frac{1}{z}g_{1}(x)\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})\right]}{\sum_{a}e^{2}_{a}f^{a}_{1}(x)D^{a}_{1,oo}(z,M_{h}^{2})}.$
(19)

Following the convention used by COMPASS in Ref. S.Sirtl , the depolarization
factors are not included in the numerator and the denominator.

##  III Calculation of the DiFF $\widetilde{D}^{\sphericalangle}_{ot}$ in the
spectator model

The twist-3 DiFF $\widetilde{D}^{\sphericalangle}$ originates from the quark-
gluon-quark correlation Bacchetta:2003vn ,

$\displaystyle\widetilde{\Delta}_{A}^{\alpha}(k,P_{h},R)$ $\displaystyle=$
$\displaystyle\frac{1}{2z}\sum_{X}\int\frac{d\xi^{+}d^{2}\xi_{T}}{(2\pi)^{3}}e^{ik\cdot\xi}\langle
0|\int_{\pm\infty^{+}}^{\xi^{+}}d\eta^{+}{\cal{U}}^{\xi_{T}}_{(\infty^{+},\xi^{+})}$
(21) $\displaystyle\times
gF_{\bot}^{-\alpha}{\cal{U}}^{\xi_{T}}_{(\eta^{+},\xi^{+})}\psi(\xi)|P_{h},R;X\rangle\langle
P_{h},R;X|\bar{\psi}(0){\cal{U}}^{0_{T}}_{(0^{+},\infty^{+})}{\cal{U}}^{\infty^{+}}_{(0_{T},\xi_{T})}|0\rangle\mid_{\eta^{+}=\xi^{+}=0,\eta_{T}=\xi_{T}}.$

Here $F_{\bot}^{-\alpha}$ is the field strength tensor of the gluon.
Introducing the covariant derivative
$iD^{\mu}(\xi)=i\partial_{\mu}+gA_{\mu}(\xi)$, we can recover also the
relation

$\displaystyle\widetilde{\Delta}_{A}^{\alpha}(k,P_{h},R)$ $\displaystyle=$
$\displaystyle{\Delta}_{D}^{\alpha}(k,P_{h},R)-{\Delta}_{\partial}^{\alpha}(k,P_{h},R)\,,$
(22)

Where

$\displaystyle{\Delta}_{D}^{\alpha}(k,P_{h},R)$ $\displaystyle=$
$\displaystyle z^{2}\sum_{X}\int\frac{d\xi^{+}}{2\pi}e^{ik\cdot\xi}\langle
0|{\cal{U}}^{+}_{[0,\xi]}\psi(\xi)iD^{\alpha}(\xi)|P_{h},R;X\rangle\langle
P_{h},R;X|\bar{\psi}(0)|0\rangle|_{\xi^{-}=\xi_{T}=0}\,,$ (23)
$\displaystyle{\Delta}_{\partial}^{\alpha}(k,P_{h},R)$ $\displaystyle=$
$\displaystyle
z^{2}k_{T}^{\alpha}\sum_{X}\int\frac{d\xi^{+}}{2\pi}e^{ik\cdot\xi}\langle
0|{\cal{U}}^{+}_{[0,\xi]}\psi(\xi)|P_{h},R;X\rangle\langle
P_{h},R;X|\bar{\psi}(0)|0\rangle|_{\xi^{-}=\xi_{T}=0}\,.$ (24)

After integrating out ${\vec{k}_{T}}$, we get

$\displaystyle\widetilde{\Delta}_{A}^{\alpha}(z,\cos\theta,M_{h}^{2},\phi_{R})=\frac{z^{2}|\vec{R}|}{8M_{h}}\int
d^{2}\vec{k}_{T}\widetilde{\Delta}_{A}^{\alpha}(k,P_{h},R).$ (25)

Then, the DiFF $\widetilde{D}^{\sphericalangle}$ can be obtained by the trace

$\displaystyle\frac{R_{T}^{\alpha}}{z}\widetilde{D}^{\sphericalangle}(z,\cos\theta,M_{h}^{2})=4\pi\textrm{Tr}[\widetilde{\Delta}_{A}^{\alpha}(z,\cos\theta,M_{h}^{2},\phi_{R})\gamma^{-}].$
(26)

Following the approach developed in Ref. Bacchetta:2006un , we will work in
the framework of a spectator model for the fragmentation process
$q\rightarrow\pi^{+}\pi^{-}X$. Here, the sum over all possible intermediate
states $X$ is replaced by an effective on-shell state-the spectator, whose
quantum numbers are the same as the initial quark and whose mass is one of the
parameters of the model. The twist-2 DiFFs $D_{1,oo}$ and
${H}^{\sphericalangle}_{1,ot}$ have been studied in Ref. Bacchetta:2006un
using the spectator model. The model was also extended to calculate the
twist-3 DiFF $\widetilde{G}^{\sphericalangle}$ in Ref. Yang:2019aan . In the
following, the calculation of the unknown DiFF
$\widetilde{D}^{\sphericalangle}_{ot}$ in the same model will be described in
details.

The corresponding diagram to calculate the quark-gluon-quark correlation for
DiFF is shown in Fig. 2, in which the left hand side corresponds to the quark-
hadron vertex $\langle P_{h},R;X|\bar{\psi}(0)|0\rangle$, and the right hand
side corresponds to the vertex $\langle
0|igF_{\perp}^{-\alpha}(\eta^{+})\psi(\xi^{+})|P_{h},R;X\rangle$ which
contains gluon rescattering.

Figure 2: Diagrammatic representation of the correlation function
$\widetilde{\Delta}_{A}^{\alpha}$ in the spectator model.

As given in details in Ref. Bacchetta:2006un , there are a few prominent
channels contributing to the $q\rightarrow\pi^{+}\pi^{-}X$ process:

1\. $q\rightarrow\pi^{+}\pi^{-}X_{1}$;

2\. $q\rightarrow\rho X_{2}\rightarrow\pi^{+}\pi^{-}X_{2}$;

3\. $q\rightarrow\omega X_{3}\rightarrow\pi^{+}\pi^{-}X_{3}$;

4\. $q\rightarrow\omega X^{\prime}_{4}\rightarrow\pi^{+}\pi^{-}X_{4}$ with
$X_{4}=\pi^{0}X^{\prime}_{4}$;

5\. $q\rightarrow\eta X^{\prime}_{5}\rightarrow\pi^{+}\pi^{-}X_{5}$ with
$X_{5}=\chi X^{\prime}_{5}$;

6\. $q\rightarrow K^{0}X_{6}\rightarrow\pi^{+}\pi^{-}X_{6}$.

In the first process, the quark fragments into an “incoherent”
$\pi^{+}\pi^{-}$ pair that are called as “background”; while in the other five
processes the $\pi^{+}pi^{-}$ pair are produced through the decays of the
intermediate resonances $\rho$, $\omega$, $\eta$ and $K^{0}$, responsible for
the peaks at $M_{h}\sim$ 770MeV, 782MeV, 500MeV, and 498MeV, respectively.

As stated in Ref. Bacchetta:2006un , different channels could produce
spectators with different masses. For simplicity, here we consider just a
single spectator for all channels. The spectator mass is denoted by $M_{s}$.
The choice of using the same spectator for all channels implies in particular
that the fragmentation amplitudes of all channels can interfere with each
other maximally. In reality, only a fraction of the total events ends up in
the same spectator and can thus produce interference effects.

Using the above setup, we can write down the quark-gluon-quark correlator for
dihardon fragmentation produced by $s$-wave and $p$-wave

$\displaystyle\widetilde{\Delta}_{A}^{\alpha}(k,P_{h},R)=i\frac{C_{F}\alpha_{s}}{2(2\pi)^{2}(1-z)P_{h}^{-}}\frac{1}{k^{2}-m^{2}}\int\frac{d^{4}l}{(2\pi)^{4}}(l^{-}g_{T}^{\alpha\mu}-l^{\alpha}_{T}g^{-\mu})$
(27)
$\displaystyle\frac{(k\\!\\!\\!/-l\\!\\!\\!/+m)(F^{s\star}e^{-\frac{k^{2}}{\Lambda_{s}^{2}}}+F^{p\star}e^{-\frac{k^{2}}{\Lambda_{p}^{2}}}R\\!\\!\\!/)(k\\!\\!\\!/-P\\!\\!\\!/_{h}-l\\!\\!\\!/+m_{s})\gamma_{\mu}(k\\!\\!\\!/-P\\!\\!\\!/_{h}+m_{s})(F^{s}e^{-\frac{k^{2}}{\Lambda_{s}^{2}}}+F^{p}e^{-\frac{k^{2}}{\Lambda_{p}^{2}}}R\\!\\!\\!/)(k\\!\\!\\!/+m)}{(-l^{-}\pm
i\epsilon)((k-l)^{2}-m^{2}-i\epsilon)((k-P_{h}-l)^{2}-m_{s}^{2}-i\epsilon)(l^{2}-i\epsilon)}.$
(28)

Here, $m$ and $m_{s}$ are the mass of the quark and the spectator,
respectively, and $\Lambda_{s}$ and $\Lambda_{p}$ are the cutoffs for the
quark momentum with the form Bacchetta:2006un

$\displaystyle\Lambda_{s,p}=\alpha_{s,p}z^{\beta_{s,p}}(1-z)^{\gamma_{s,p}},$
(29)

The exponential form factors are introduced to suppress the contributions from
higher quark virtuality Gamberg:2003eg . The open circle in Fig. 2 denote The
factor $(l^{-}g_{T}^{\alpha\mu}-l_{T}^{\alpha}g^{-\mu})$ coming form the
Feynman rule for the gluon field strength tensor in the operator definition of
the correlator 21. In addition, $F^{s}$ and $F^{p}$ are the quark-dihadron-
spectator vertices associated with the $s$-wave and $p$-wave contributions,
particular, the $p$-wave vertices contain the real part and image part, which
are both contribution to the twist-3 fragmentation function
$\widetilde{D}^{\sphericalangle}$, respectively, and in Ref. Bacchetta:2006un
they are parameterized as:

$\displaystyle F^{s}$ $\displaystyle=f_{s}\,,$ $\displaystyle F^{p}$
$\displaystyle=f_{\rho}\frac{(M_{h}^{2}-M_{\rho}^{2})-i\Gamma_{\rho}M_{\rho}}{(M_{h}^{2}-M_{\rho}^{2})+\Gamma_{\rho}^{2}M_{\rho}^{2}}+f_{\omega}\frac{(M_{h}^{2}-M_{\omega}^{2})-i\Gamma_{\omega}M_{\omega}}{(M_{h}^{2}-M_{\rho}^{2})+\Gamma_{\omega}^{2}M_{\omega}^{2}}$
$\displaystyle-
if^{\prime}_{\omega}\frac{\sqrt{\hat{\lambda}(M_{\omega}^{2},M_{h}^{2},m_{\pi}^{2})}\Theta(M_{\omega}-m_{\pi}-M_{h})}{4\pi\Gamma_{\omega}^{2}[4M_{\omega}^{2}m_{\pi}^{2}+\hat{\lambda}(M_{\omega}^{2},M_{h}^{2},m_{\pi}^{2})]^{\frac{1}{4}}}\,.$
(30)

Here,
$\hat{\lambda}(M_{\omega}^{2},M_{h}^{2},m_{\pi}^{2})=(M_{\omega}^{2}-(M_{h}+m_{\pi})^{2})(M_{\omega}^{2}-(M_{h}-m_{\pi})^{2})$
and $\Theta$ denotes the unit step function. The first two terms of $F^{p}$
can be identified with the contributions of the $\rho$ and the $\omega$
resonances decaying into two pions. The masses and widths of the two
resonances are taken from the PDG Eidelman : $M_{\rho}=0.776\mathrm{GeV}$,
$\Gamma_{\rho}=0.150\mathrm{GeV}$, $M_{\omega}=0.783\mathrm{GeV}$,
$\Gamma_{\omega}=0.008\mathrm{GeV}$.

To identify the contributions from the correlator to
$\widetilde{D}^{\sphericalangle}_{ot}$, we expand (27) in the following way,

$\displaystyle\widetilde{\Delta}_{A}^{\alpha}(z,\cos\theta,M_{h}^{2},\phi_{R})=i\frac{C_{F}\alpha_{s}z^{2}|\vec{R}|}{16(2\pi)^{5}(1-z)M_{h}P_{h}^{-}}\int
d|\vec{k}_{T}|^{2}\int
d^{4}l\frac{l^{-}g_{T}^{\alpha\mu}-l^{\alpha}_{T}g^{-\mu}}{k^{2}-m^{2}}$ (31)
$\displaystyle\bigg{[}+|F^{s}|e^{-\frac{2k^{2}}{\Lambda_{s}^{2}}}\frac{(k\\!\\!\\!/-l\\!\\!\\!/+m)(k\\!\\!\\!/-P\\!\\!\\!/_{h}-l\\!\\!\\!/+m_{s})\gamma_{\mu}(k\\!\\!\\!/-P\\!\\!\\!/_{h}+m_{s})(k\\!\\!\\!/+m)}{(-l^{-}\pm
i\epsilon)((k-l)^{2}-m^{2}-i\epsilon)((k-P_{h}-l)^{2}-m_{s}^{2}-i\epsilon)(l^{2}-i\epsilon)}$
(32)
$\displaystyle+|F^{p}|e^{-\frac{2k^{2}}{\Lambda_{p}^{2}}}\frac{(k\\!\\!\\!/-l\\!\\!\\!/+m)R\\!\\!\\!/(k\\!\\!\\!/-P\\!\\!\\!/_{h}-l\\!\\!\\!/+m_{s})\gamma_{\mu}(k\\!\\!\\!/-P\\!\\!\\!/_{h}+m_{s})R\\!\\!\\!/(k\\!\\!\\!/+m)}{(-l^{-}\pm
i\epsilon)((k-l)^{2}-m^{2}-i\epsilon)((k-P_{h}-l)^{2}-m_{s}^{2}-i\epsilon)(l^{2}-i\epsilon)}$
(33)
$\displaystyle+(F^{s\star}F^{p})e^{-\frac{2k^{2}}{\Lambda_{sp}^{2}}}\frac{(k\\!\\!\\!/-l\\!\\!\\!/+m)(k\\!\\!\\!/-P\\!\\!\\!/_{h}-l\\!\\!\\!/+m_{s})\gamma_{\mu}(k\\!\\!\\!/-P\\!\\!\\!/_{h}+m_{s})R\\!\\!\\!/(k\\!\\!\\!/+m)}{(-l^{-}\pm
i\epsilon)((k-l)^{2}-m^{2}-i\epsilon)((k-P_{h}-l)^{2}-m_{s}^{2}-i\epsilon)(l^{2}-i\epsilon)}$
(34)
$\displaystyle+(F^{s}F^{p\star})e^{-\frac{2k^{2}}{\Lambda_{sp}^{2}}}\frac{(k\\!\\!\\!/-l\\!\\!\\!/+m)R\\!\\!\\!/(k\\!\\!\\!/-P\\!\\!\\!/_{h}-l\\!\\!\\!/+m_{s})\gamma_{\mu}(k\\!\\!\\!/-P\\!\\!\\!/_{h}+m_{s})(k\\!\\!\\!/+m)}{(-l^{-}\pm
i\epsilon)((k-l)^{2}-m^{2}-i\epsilon)((k-P_{h}-l)^{2}-m_{s}^{2}-i\epsilon)(l^{2}-i\epsilon)}\bigg{]}\,.$
(35)

where $2/\Lambda^{2}_{sp}=1/\Lambda^{2}_{s}+1/\Lambda^{2}_{p}$.

In the spectator model, the factor $k^{2}$ in Eq. (31), compatible with the
on-shell condition of the spectator, is given by

$\displaystyle k^{2}$ $\displaystyle=$
$\displaystyle\frac{z}{1-z}|\vec{k}_{T}|^{2}+\frac{M_{s}^{2}}{1-z}+\frac{M_{h}^{2}}{z}\,.$
(36)

In Eq. (31), the first term corresponds to the pure $s$-wave contribution, the
second term corresponds to the pure $p$-wave contribution; whereas the third
and fourth term correspond to the interference of $s$ and $p$-wave
contributions. According to the partial-wave analysis for
$\widetilde{D}^{\sphericalangle}$, one can find that only the third and fourth
terms contribute to $\widetilde{D}^{\sphericalangle}_{ot}$. Integrating over
the internal momentum $l$ yields the final expression for
$\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})$:

$\displaystyle\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})$
$\displaystyle=-\frac{\alpha_{s}C_{F}z^{2}|\vec{R}|}{4(2\pi)^{4}(1-z)M_{h}}\int
d|\vec{k}_{T}|^{2}e^{-\frac{2k^{2}}{\Lambda_{sp}^{2}}}\frac{1}{k^{2}-m^{2}}\bigg{\\{}\textrm{Re}(F^{s*}F^{p})C$
$\displaystyle+\textrm{Im}(F^{s*}F^{p})m_{s}\big{[}(k^{2}-m^{2})(A+zB)-2(Ak^{2}+BP_{h}\cdot
k)\big{]}\bigg{\\}}\,,$ (37)

Here, the coefficient $C$ corresponds to real part of the integration over $l$
and has the expression

$\displaystyle
C=m\int_{0}^{1}dx\int_{0}^{1-x}dy\frac{[x+(1-z)y-2](k^{2}-2k\cdot
P_{h}-m_{s}^{2}+m_{h}^{2})+(x+2y)k\cdot
P_{h}-(x+y)k^{2}-ym_{h}^{2}}{x(1-x)k^{2}+2k\cdot(k-P_{h})xy+m^{2}x+m_{s}^{2}y+y(y-1)(k-P_{h})^{2}}\,,$
(38)

Note that it is proportional to the fragmenting quark mass $m$. The factor
${A}$ and ${B}$ appear in Eq. (37) have the forms

$\displaystyle A$
$\displaystyle={I_{1}\over\lambda(m_{h},m_{s})}\left(2k^{2}\left(k^{2}-m_{s}^{2}-m_{h}^{2}\right){I_{2}\over\pi}+\left(k^{2}+m_{h}^{2}-m_{s}^{2}\right)\right),$
(39) $\displaystyle B$
$\displaystyle=-{2k^{2}\over\lambda(m_{h},m_{s})}I_{1}\left(1+{k^{2}+m_{s}^{2}-m_{h}^{2}\over\pi}I_{2}\right)\,.$
(40)

And the functions $I_{i}$ have the forms Amrath:2005gv

$\displaystyle I_{1}$ $\displaystyle=\int
d^{4}l\delta(l^{2})\delta((k-l)^{2}-m^{2})={\pi\over
2k^{2}}\left(k^{2}-m^{2}\right)\,,$ (41) $\displaystyle I_{2}$
$\displaystyle=\int
d^{4}l{\delta(l^{2})\delta((k-l)^{2}-m^{2})\over(k-P_{h}-l)^{2}-m_{s}^{2}}={\pi\over
2\sqrt{\lambda(m_{h},m_{s})}}\ln\left(1-{2\sqrt{\lambda(m_{h},m_{s})}\over
k^{2}-m_{h}^{2}+m_{s}^{2}+\sqrt{\lambda(m_{h},m_{s})}}\right)\,.$ (42)

with
$\lambda(m_{h},m_{s})=(k^{2}-(m_{h}+m_{s})^{2})(k^{2}-(m_{h}-m_{s})^{2})$.

##  IV Numerical Estimate for $\widetilde{D}^{\sphericalangle}_{ot}$ and the
double spin asymmetry

In order to obtain the numerical results of the fragmentation function
$\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})$, we need to know the
values for the model parameters $m_{s}$, $\alpha_{s,p}$, $\beta_{s,p}$ and
$\gamma_{s,p}$. For these we adopt them from Ref. Bacchetta:2006un :

$\displaystyle\alpha_{s}$ $\displaystyle=2.60\pm
0.05\mathrm{GeV}^{2}\,,\qquad\beta_{s}=-0.751\pm
0.008\,,\qquad\gamma_{s}=-0.193\pm 0.004\,,$ $\displaystyle\alpha_{p}$
$\displaystyle=7.07\pm 0.11\mathrm{GeV}^{2}\,,\qquad\beta_{p}=-0.038\pm
0.003\,,\qquad\gamma_{p}=-0.085\pm 0.004\,,$ $\displaystyle f_{s}$
$\displaystyle=1197\pm 2\mathrm{GeV}^{-1}\,,\qquad f_{\rho}=93.5\pm
1.6\,,\qquad f_{\omega}=0.63\pm 0.03\,,$ $\displaystyle f^{\prime}_{\omega}$
$\displaystyle=75.2\pm 1.2\,,\qquad M_{s}=2.97\pm 0.04M_{h}\,.$ (43)

For the quark mass $m$, we will adopt two different choices for comparison.
The first one is $m=0$ GeV following the adoption in Ref. Bacchetta:2006un .
In this choice the $C$ term (Eq. (38)) vanishes since it is proportional to
$m$. Therefore, in this particular case only the term containing
$\textrm{IM}(F^{s*}F^{p})$ in Eq. (37) contributes to
$\widetilde{D}^{\sphericalangle}_{ot}$ numerically. In the second choice we
adopt $m=0.3$ GeV, which is consistent with the value for the quark massl
chosen in Refs. Bacchetta:2008af  and Bacchetta:2007wc .

The results for $\widetilde{D}^{\sphericalangle}_{ot}$ divided by the
unpolarized DiFF $D_{1,oo}$ are plotted in Fig. 3. The left panel depicts the
$z$-dependent of the ratio ($M_{h}$ is integrated out in the region
$0.3<M_{h}<1.6$), and the right panel depicts the $M_{h}$-dependence of the
ratio ($z$ is integrated out over the region $0.2<z<0.9$). The dashed line
corresponds to the result for $m=0$ GeV. The dashed-dotted line corresponds to
the result from the Im$(F^{s*}F^{p})$ term for $m=0.3$ GeV, while the dotted
line corresponds to the nonzero result from the $C$ term in Eq. (38) for
$m=0.3$ GeV. The solid line is the total result for $m=0.3$ GeV, corresponding
to the sum of the dashed-dotted and dotted lines. We find that
${\widetilde{D}^{\sphericalangle}_{ot}}/{D_{1,oo}}$ is positive in the entire
$z$ and $M_{h}$ region when $M_{h}$ or $z$ is integrated out, respectively. We
note that the $\cos R_{\phi}$ asymmetry in the unpolarized SIDIS is
proportional to the ratio $\widetilde{D}^{\sphericalangle}/D_{1}$ as suggested
by Eq. (44) in Ref. Bacchetta:2003vn , therefore, the curve in Fig. 3 can be
also viewed as an approximate result for the $\cos\phi_{R}$ asymmetry in the
unpolarized SIDIS. We also observe that two difference choices on the quark
mass $m$ lead to a $10\%$ difference on the ratio
${\widetilde{D}^{\sphericalangle}_{ot}}/{D_{1,oo}}$

Figure 3: The twist-3 D0iFF $\widetilde{D}^{\sphericalangle}_{ot}$ as the
functions of $z$ (left panel) and $M_{h}$ (right panel) in the spectator
model, normalized by the unplarized DiFF $D_{1,oo}$. The dashed line
corresponds to the result for $m=0$ GeV. The dashed-dotted line corresponds to
the result from the Im$(F^{s*}F^{p})$ term for $m=0.3$ GeV, while the dotted
line corresponds to the result from the $C$ term in Eq. (38) for $m=0.3$ GeV.
The solid line is the total result for $m=0.3$ GeV, corresponding to the sum
of the dashed-dotted and dotted lines.

Using Eq. (19), one can express the double longitudinally polarized asymmetry
$A^{\cos\phi_{R}}_{LL}$ in SIDIS as functions of $x$, $z$ and $M_{h}$ as
follow

$\displaystyle A^{\cos\phi_{R}}_{LL}(x)=-\frac{\int dz\int
dM_{h}2M_{h}\frac{|\vec{R}|}{Q}\frac{1}{z}(4g_{1}^{u}(x)+g_{1}^{d}(x))\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})}{\int
dz\int dM_{h}2M_{h}(4f_{1}^{u}(x)+f_{1}^{d}(x))D_{1,oo}(z,M_{h}^{2})},$ (44)
$\displaystyle A^{\cos\phi_{R}}_{LL}(z)=-\frac{\int dx\int
dM_{h}2M_{h}\frac{|\vec{R}|}{Q}\frac{1}{z}(4g_{1}^{u}(x)+g_{1}^{d}(x))\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})}{\int
dx\int dM_{h}2M_{h}(4f_{1}^{u}(x)+f_{1}^{d}(x))D_{1,oo}(z,M_{h}^{2})},$ (45)

$\displaystyle A^{\cos\phi_{R}}_{LL}(M_{h})=-\frac{\int dx\int
dz\frac{|\vec{R}|}{Q}\frac{1}{z}(4g_{1}^{u}(x)+g_{1}^{d}(x))\widetilde{D}^{\sphericalangle}_{ot}(z,M_{h}^{2})}{\int
dx\int dz(4f_{1}^{u}(x)+f_{1}^{d}(x))D_{1,oo}(z,M_{h}^{2})}.$ (46)

For the unpolarized distribution $f_{1}$ and the helicity distribution
$g_{1}$, we also apply a spectator model result from Ref. Bacchetta:2008af
for consistency. As the scale dependence of the DiFF
$\widetilde{D}^{\sphericalangle}$ still remains unknown, we assume that
$\widetilde{D}^{\sphericalangle}$ follows the same evolution as that of the
DiFF $D_{1}$ Ceccopieri:2007ip ; deFlorian:2003cg  in leading order. For the
PDFs $f_{1}$ and $g_{1}$, we adopt the leading-order DGLAP evolution. We admit
that the scale dependence of $\widetilde{D}^{\sphericalangle}$ might be more
complicated than that of $D_{1}$ because $\widetilde{D}^{\sphericalangle}$ is
a twist-3 object. Since here the asymmetry is the ratio between
$g_{1}\otimes\widetilde{D}^{\sphericalangle}$ and $f_{1}\otimes{D}_{1}$, we
supposed the evolution effect will not influence the results qualitatively.

The COMPASS Collaboration S.Sirtl  is measuring the above $A_{LL}^{\cos\phi}$
asymmetry using a 160 GeV or a 190 GeV longitudinally polarized muon beam on a
longitudinally polarized nucleon target. Using the numerical result for
$\widetilde{D}^{\sphericalangle}_{ot}$, we estimate the asymmetry
$A_{LL}^{\cos\phi_{R}}$ at the kinematical region of COMPASS

$\displaystyle 0.003<x<0.4,\quad 0.1<y<0.9,\quad 0.2<z<0.9,$ $\displaystyle
0.3\mathrm{GeV}<M_{h}<1.6\mathrm{GeV},\quad 1\mathrm{GeV}<Q^{2},\quad
W>5\mathrm{GeV}.$ (47)

Here $W$ is invariant mass of the virtual photon-nucleon system. The
kinematical cuts in Eq. (47) are implemented in the calculation numerically by
imposing those constraints in the computation code.

In the left, central, and right panels of Fig. 4, we plot the $x$, $z$ and
$M_{h}$ dependent $\cos\phi_{R}$ asymmetry at the COMPASS kinematics. The
solid and dashed lines correspond to the results for $m=0.3$ GeV ad $m=0$ GeV,
respectively. Again, the difference between this two results is about $10\%$.
We find that the asymmetry is negative due to the minus sign in Eq. (19). The
asymmetry roughly decreases with increasing $x$. The estimated magnitude of
the asymmetry at COMPASS is about 0.01, which is smaller than the
$\sin\phi_{R}$ asymmetry in single longitudinally polarized SIDIS S.Sirtl ;
Yang:2019aan .

Figure 4: The $\cos\phi_{R}$ azimuthal asymmetry in dihadron production off
the longitudinally polarized proton as functions of $x$ (left panel), $z$
(central panel) and $M_{h}$ (right panel) at COMPASS. The dashed and solid
lines corresponds to the results from $m=0$ GeV and $m=0.3$ GeV.

For comparison, we also predict the asymmetry $A_{LL}^{\cos\phi_{R}}$ at JLab
12GeV using the following kinematical configuration

$\displaystyle\quad 0.072<x<0.532,\quad 0.2<y<0.95,\quad 0.2<z<0.8,\quad
0.5\textrm{GeV}<M_{h}<1.2\textrm{GeV},$ $\displaystyle\quad
E_{e}=12\mathrm{GeV},\quad W>4\mathrm{GeV}^{2},\quad
1<Q^{2}<6.3\mathrm{GeV}^{2}.$ (48)

The corresponding results are shown in Fig. 5. In general, The asymmetry at
JLab 12GeV is found to be larger than that at COMPASS because the
$A_{LL}^{\cos\phi_{R}}$ asymmetry is a twist-3 effect and the c.m. energy at
JLab is smaller.

Figure 5: The $\cos\phi_{R}$ azimuthal asymmetry in dihadron production off
the longitudinally polarized proton as functions of $x$ (left panel), $z$
(central panel) and $M_{h}$ (right panel) at JLab 12GeV. The dashed and solid
lines corresponds to the results from $m=0$ GeV and $m=0.3$ GeV.

We also present the $\cos\phi_{R}$ asymmetry in double polarized SIDIS at the
kinematical configuration of a future EIC facility Matevosyan:2015gwa :

$\displaystyle\sqrt{s}=45\mathrm{GeV},\quad 0.001<x<0.4,\quad
0.01<y<0.95,\quad 0.2<z<0.8,$ $\displaystyle
0.3\textrm{GeV}<M_{h}<1.6\textrm{GeV},\quad Q^{2}>1\mathrm{GeV}^{2},\quad
W>5\mathrm{GeV}.$ (49)

The asymmetry vs $x$, $z$ and $M_{h}$ are plotted in the left, central and
right panels of Fig. 6. We find that the overall tendency of the asymmetry at
the EIC is similar to that at COMPASS. Although the size of the asymmetry is
smaller due to the higher-twist nature of the asymmetry, it is still
measurable at the kinematics of EIC.

Figure 6: The $\cos\phi_{R}$ azimuthal asymmetry in dihadron production off
the longitudinally polarized proton as functions of $x$ (left panel), $z$
(central panel) and $M_{h}$ (right panel) at an EIC. The dashed and solid
lines corresponds to the results from $m=0$ GeV and $m=0.3$ GeV.

Finally, we note that a recent measurement on the single-transverse spin
asymmetry of dihadron production in SIDIS by COMPASS Adolph:2014fjw  shows
that the asymmetry is similar to the Collins asymmetry of single hadron
production in SIDIS. This hints that the spin-dependent DiFFs could be
generated by the single-hadron spin-dependent fragmentation functions. The
closed relation between the Collins fragmentation function $H_{1}^{\perp}$ and
the spin-dependent DiFF $H_{1}^{\sphericalangle}$ in large $M_{h}$ region has
also been suggested in Ref. Zhou:2011ba , which showed that both the leading
twist unpolarized dihadron fragmentation function and the leading twist
interference FF (IFF) can be perturbatively computed in terms of the normal
single hadron fragmentation function when the invariant mass is sufficiently
large. Therefore, there is the possibility that the similar analysis can be
applied to the twist-3 IFF case as well. On the other hand, a calculation of
SSA in dihadron production based on the spectator model result
Bacchetta:2006un  for $H_{1}^{\sphericalangle}$ can also well describe the
COMPASS data (see solid blue lines in Fig. 5 of Ref. Adolph:2014fjw ) for the
$x$ and $z$ dependent shape, as well as for the $M_{h}$ dependent shape in the
$\rho$ meson region. This may suggest that the resonance interference
mechanism applied in our paper can also responsible for the spin asymmetry in
dihadron production. Due to the lack of further data, at the moment it is hard
to say which mechanism should be preferred. Further theoretical and
experimental studies are needed in order to discriminate different mechanisms
for the spin-dependent DiFFs.

##  V Conclusion

In this work, we studied the origin of the $\cos\phi_{R}$ asymmetry of hadron
pair production in double polarized SIDIS:
$l^{\rightarrow}+p^{\rightarrow}\rightarrow h_{1}+h_{2}+X$ The asymmetry can
originate from the coupling of the twist-3 DiFF
$\widetilde{D}^{\sphericalangle}$ and the helicity distribution $g_{1}(x)$.
Another potential contribution, the coupling
$e_{L}(x)\,H_{1}^{\sphericalangle}$, will not give rise to the asymmetry
because of the time-reversal invariant constraint $\int
d^{2}p_{T}\,e_{L}(x,p^{2}_{T})=0$. We applied a spectator model for the quark-
gluon-quark correlator and calculated the $s$-wave and $p$-wave interference
DiFF $\widetilde{D}^{\sphericalangle}_{ot}$, the leading term of
$\widetilde{D}^{\sphericalangle}$ in the partial wave expansion. Using the
numerical results for $\widetilde{D}^{\sphericalangle}_{ot}$, we estimated the
$\cos\phi_{R}$ asymmetry as functions of $x$, $z$ and $M_{h}$ at the
kinematics of COMPASS, JLab 12GeV, and EIC. We found that the asymmetry at
COMPASS and JLab12 is about 1-2 percent and may be measurable. Therefore, the
measurement of the $\cos\phi_{R}$ asymmetry at COMPASS and JLab12 may provide
unambiguous information on $\widetilde{D}^{\sphericalangle}$.

##  VI Acknowledgements

This work is partially supported by the NSFC (China) grant 11575043. W. Y and
H. L contributed equally to this work and should be considered as co-first
authors.

## References

  * (1) P. J. Mulders and R. D. Tangerman, Nucl. Phys. B 461, 197 (1996) Erratum: [Nucl. Phys. B 484, 538 (1997)]
  * (2) A. Kotzinian, Nucl. Phys. B 441, 234 (1995)
  * (3) D. Boer and P. J. Mulders, Phys. Rev. D 57, 5780 (1998) [hep-ph/9711485].
  * (4) K. Goeke, A. Metz and M. Schlegel, Phys. Lett. B 618, 90 (2005) doi:10.1016/j.physletb.2005.05.037 [hep-ph/0504130].
  * (5) A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders and M. Schlegel, JHEP 0702, 093 (2007)
  * (6) X. d. Ji, J. p. Ma and F. Yuan, Phys. Rev. D 71, 034005 (2005)
  * (7) A. V. Efremov, K. Goeke and P. Schweitzer, Phys. Rev. D 67, 114014 (2003) [hep-ph/0208124].
  * (8) L. P. Gamberg, D. S. Hwang and K. A. Oganessyan, Phys. Lett. B 584, 276 (2004) [hep-ph/0311221].
  * (9) H. Avakian, A. V. Efremov, K. Goeke, A. Metz, P. Schweitzer and T. Teckentrup, Phys. Rev. D 77, 014023 (2008) [arXiv:0709.3253 [hep-ph]].
  * (10) W. Mao and Z. Lu, Phys. Rev. D 87, no. 1, 014012 (2013) [arXiv:1210.4790 [hep-ph]].
  * (11) W. Mao, Z. Lu and B. Q. Ma, Phys. Rev. D 90, no. 1, 014048 (2014) [arXiv:1405.3876 [hep-ph]].
  * (12) W. Mao, Z. Lu, B. Q. Ma and I. Schmidt, Phys. Rev. D 91, no. 3, 034029 (2015) [arXiv:1412.7390 [hep-ph]].
  * (13) Y. k. Song, J. h. Gao, Z. t. Liang and X. N. Wang, Phys. Rev. D 83, 054010 (2011) [arXiv:1012.4179 [hep-ph]].
  * (14) S. y. Wei, Y. k. Song, K. b. Chen and Z. t. Liang, Phys. Rev. D 95, no. 7, 074017 (2017) [arXiv:1611.08688 [hep-ph]].
  * (15) Y. Yang, Z. Lu and I. Schmidt, Phys. Lett. B 761, 333 (2016).
  * (16) C. Lorcé, B. Pasquini and P. Schweitzer, Eur. Phys. J. C 76, no. 7, 415 (2016) [arXiv:1605.00815 [hep-ph]].
  * (17) S. Bastami et al., arXiv:1807.10606 [hep-ph].
  * (18) A. Bacchetta, D. Boer, M. Diehl and P. J. Mulders, JHEP 0808 (2008) 023 [arXiv:0803.0227 [hep-ph]].
  * (19) L. P. Gamberg, D. S. Hwang, A. Metz and M. Schlegel, Phys. Lett. B 639, 508 (2006) [hep-ph/0604022].
  * (20) A. Bacchetta and M. Radici, Phys. Rev. D 69, 074026 (2004)
  * (21) A. Airapetian et al. [HERMES Collaboration], Phys. Lett. B 648, 164 (2007) [hep-ex/0612059].
  * (22) A. Airapetian et al. [HERMES Collaboration], Phys. Lett. B 622, 14 (2005) [hep-ex/0505042].
  * (23) M. G. Alekseev et al. [COMPASS Collaboration], Eur. Phys. J. C 70, 39 (2010) [arXiv:1007.1562 [hep-ex]].
  * (24) M. Aghasyan et al., Phys. Lett. B 704, 397 (2011) [arXiv:1106.2293 [hep-ex]].
  * (25) W. Gohn et al. [CLAS Collaboration], Phys. Rev. D 89, no. 7, 072011 (2014) [arXiv:1402.4097 [hep-ex]].
  * (26) A. Airapetian et al. [HERMES Collaboration], arXiv:1903.08544 [hep-ex].
  * (27) S. Sirtl, in 22nd International Symposium on Spin Physics (SPIN 2016) Urbana, IL, USA (unpublished) arXiv: 1702.07317, http://www.indiana.edu/ spin2016/.
  * (28) A. Bacchetta, A. Courtoy and M. Radici, Phys. Rev. Lett 107, 012001 (2011) [arXiv:1104.3855 [hep-ph]].
  * (29) A. Bacchetta, A. Courtoy and M. Radici, J. High Energy Phys 1303, (2013) 119 [arXiv:1212.3568 [hep-ph]].
  * (30) M. Radici, A. Courtoy, A. Bacchetta and M. Guagnelli, J. High Energy Phys 1505, (2015) 123 [arXiv:1503.03495 [hep-ph]].
  * (31) M. Radici and A. Bacchetta, Phys. Rev. Lett 120, 192001 (2018)
  * (32) A. Airapetian et al. [HERMES Collaboration], JHEP 0806, 017 (2008) [arXiv:0803.2367 [hep-ex]].
  * (33) C. Adolph et al. [COMPASS Collaboration], Phys. Lett. B 713, 10 (2012) [arXiv:1202.6150 [hep-ex]].
  * (34) C. Adolph et al. [COMPASS Collaboration], Phys. Lett. B 736, 124 (2014).
  * (35) J. C. Collins, S. F. Heppelmann and G. A. Ladinsky, Nucl. Phys. B 420, 565 (1994)
  * (36) R. L. Jaffe, X. m. Jin and J. Tang, Phys. Rev. Lett. 80, 1166 (1998)
  * (37) A. Bacchetta and M. Radici, Phys. Rev. D 67, 094002 (2003)
  * (38) M. Boglione and P. J. Mulders, Phys. Rev. D 60, 054007 (1999)
  * (39) D. Boer, Phys. Rev. D 60, 014012 (1999)
  * (40) A. Bianconi, S. Boffi, R. Jakob and M. Radici, Phys. Rev. D 62, 034008 (2000) [hep-ph/9907475].
  * (41) M. Radici, R. Jakob and A. Bianconi, Phys. Rev. D 65, 074031 (2002) [hep-ph/0110252].
  * (42) A. Courtoy, A. Bacchetta, M. Radici and A. Bianconi, Phys. Rev. D 85, 114023 (2012) [arXiv:1202.0323 [hep-ph]].
  * (43) A. Bacchetta and M. Radici, Phys. Rev. D 74, 114007 (2006)
  * (44) D. Boer, P. J. Mulders and F. Pijlman, Nucl. Phys. B 667, 201 (2003)
  * (45) W. Yang, X. Wang, Y. Yang and Z. Lu, Phys. Rev. D 99, no. 5, 054003 (2019) [arXiv:1902.07889 [hep-ph]].
  * (46) L. P. Gamberg, G. R. Goldstein and K. A. Oganessyan, Phys. Rev. D 68, 051501 (2003)
  * (47) S. Eidelman et al. (Particle Data Group), Phys. Lett. B592, 1 (2004).
  * (48) D. Amrath, A. Bacchetta and A. Metz, Phys. Rev. D 71, 114018 (2005), arXiv:hep-ph/0504124.
  * (49) A. Bacchetta, F. Conti and M. Radici, Phys. Rev. D 78, 074010 (2008) [arXiv:0807.0323 [hep-ph]].
  * (50) A. Bacchetta, L. P. Gamberg, G. R. Goldstein and A. Mukherjee, Phys. Lett. B 659, 234 (2008) [arXiv:0707.3372 [hep-ph]].
  * (51) F. A. Ceccopieri, M. Radici and A. Bacchetta, Phys. Lett. B 650, 81 (2007) [hep-ph/0703265 [HEP-PH]].
  * (52) D. de Florian and L. Vanni, Phys. Lett. B 578, 139 (2004) [hep-ph/0310196].
  * (53) H. H. Matevosyan, A. Kotzinian, E. C. Aschenauer, H. Avakian and A. W. Thomas, Phys. Rev. D 92, 054028 (2015) [arXiv:1502.02669 [hep-ph]].
  * (54) J. Zhou and A. Metz, Phys. Rev. Lett. 106, 172001 (2011) [arXiv:1101.3273 [hep-ph]].


# Lipschitz continuity results
for a class of obstacle problems

Carlo Benassi, Michele Caselli  Dipartimento di Scienze Fisiche, Informatiche
e Matematiche, via Campi 213/b, 41125 Mo-dena, Italy

###### Abstract.

We prove Lipschitz continuity results for solutions to a class of obstacle
problems under standard growth conditions of $p-$type, $p\geq 2$. The main
novelty is the use of a linearization technique going back to [28] in order to
interpret our constrained minimizer as a solution to a nonlinear elliptic
equation, with a bounded right hand side. This lead us to start a Moser
iteration scheme which provides the $L^{\infty}$ bound for the gradient. The
application of a recent higher differentiability result [24] allows us to
simplify the procedure of the identification of the Radon measure in the
linearization technique employed in [32]. To our knowdledge, this is the first
result for non-automonous functionals with standard growth conditions in the
direction of the Lipschitz regularity.

###### Key words and phrases:

Regularity of solutions, obstacle problems, a priori estimates

###### 2000 Mathematics Subject Classification:

Primary: 49N60, 35J85; Secondary: 35B65, 35B45, 49J40

The authors are partially supported GNAMPA (Gruppo Nazionale per l’Analisi
Matematica, la Probabilità e le loro Applicazioni) of INdAM (Istituto
Nazionale di Alta Matematica) and by the University of Modena and Reggio
Emilia through the project FAR2017 “Equazioni differenziali: problemi
evolutivi, variazionali ed applicazioni” (coord. Prof. S. Gatti). Part of the
research has been performed while the second author was visiting the
University of Napoli “Federico II” and the University of Napoli Parthenope.
The hospitality of these Institutions, in particular of Prof. Antonia
Passarelli di Napoli and Prof. Raffaella Giova are warmly acknowledged.

##  1\. Introduction

The aim of this paper is the study of the Lipschitz continuity of the
solutions $u\in W^{1,p}(\Omega)$ to variational obstacle problems of the form

$\min\left\\{\int_{\Omega}f(x,Dw):w\in\mathcal{K}_{\psi,\Psi}(\Omega)\right\\}.$
(1.1)

The function $\psi:\,\Omega\rightarrow[-\infty,+\infty)$, called obstacle,
belongs to the Sobolev space $W^{1,p}(\Omega)$ and the class
$\mathcal{K}_{\psi,\Psi}(\Omega)$ is defined as follows, for a given $\Psi\in
W^{1,p}(\overline{\Omega})$ providing the boundary value

$\mathcal{K}_{\psi,\Psi}(\Omega):=\left\\{w\in
W^{1,p}(\Omega):w(x)\geq\psi(x)\,\,\textnormal{a.e. in $\Omega$ and
}\,w|_{\partial\Omega}=\Psi|_{\partial\Omega}\right\\}.$ (1.2)

In the following, we shall assume that $\mathcal{K}_{\psi,\Psi}(\Omega)$ is
not empty.

The study of the regularity theory for obstacle problems is a classical topic
in Calculus of Variations and Partial Differential Equations. The first
results concerning obstacle problems date back to the sixties to the
pioneering work by G. Stampacchia [50] and G. Fichera [26]. The well known
fact that solutions to the obstacle problem cannot be of class
$\mathcal{C}^{2}$ independently of the regularity of the obstacle led to the
origin of the concept of weak solution and to the theory of variational
inequalities, after the fundamental work of J.L. Lions and G. Stampacchia
[39]. These problems can generally be solved by means of nowadays classical
methods of functional analysis; in this respect the question to establish
conditions so that weak solutions can be actually classical ones is of crucial
importance (see [8]). For more details we refer to the following monographs
[1], [10], [27], [37], [49].

The regularity of solutions to the obstacle problems is influenced by the one
of the obstacle: in the linear case, for instance, obstacle and solutions
share the same regularity ([6], [9], [37]), while in the nonlinear setting the
situation is more involved. This generated along the years an intense research
activity aimed to establish the regularity of the solution compared with the
regularity assumed for the obstacle. In particular many results concern the
Hölder continuity of solutions to the obstacle problem when the obstacle
itself is Hölder continuous: see for instance [40], [11], [17], [18], [20],
[7], [46]; see also [12], [38], [32], [3], [19], [21], [47], [2], [48], [4],
[5], [24], [25].

As far as we know, such analysis has not been carried out in the direction of
obtaining Lipschitz regularity results. The reason may be found in the
difficulty of choosing the proper test function, in order to start the
iteration process, necessary to get the desired regularity; this nowadays
classical technique ([45]) has been adapted along the years in several
situations, also in the case of non standard growth conditions, see also for
instance [41], [42], [43], [44].

In [32] Fuchs and Mingione considered constrained local minimizers of elliptic
variational integrals, assuming for the Lagrangian $f$ nearly linear growth.
In this respect, they were able to prove several regularity results; the one
we are mainly interested in is the Lipschitz continuity of these constrained
minimizers. This step has been achived through a linearization technique,
going back to [28], later refined in [16], see also [29], [30], [31].

The main idea of the linearization approach is the following: the constrained
minimizer is interpreted as a solution to an elliptic equation with a bounded
right hand side, which is obtained after the identification of a suitable
Radon measure. At this point the main obstruction is overcome, and higher
regularity for the solution to the obstacle problem can be shown.

We would like to remark that in [32] a lot of effort has been employed to
identify the Radon measure and the authors explicitly say that this procedure
could be significantly simplified if we would have a priori proved higher
differentiability for local minimizers of the obstacle problem.

However this is exactly the case for our specific situation: in a very recent
paper [24] Eleuteri and Passarelli di Napoli were able to establish the higher
differentiability of integer and fractional order of the solutions to a class
of obstacle problems (involving $p-$harmonic operators) assuming that the
gradient of the obstacle possesses an extra (integer or fractional)
differentiability property. We will make use of this result to simplify the
procedure outlined in [32] and obtain the desired result.

The paper is organized as follows: Section 2 contains some notations and the
statement of the main result; Section 3 contains the detailed construction of
the linearization procedure while the conclusion comes with Section 4 where
the a priori estimate is presented.

We would like to mention that, in order to let the linearization technique be
the main focus of our work, we assumed standard growth conditions for the
lagrangians of our integral functionals and $L^{\infty}$ dependence with
respect to the $x-$variable. In a similar way, we assume on the gradient of
the obstacle $W^{1,\infty}$ regularity. We expect that future improvements can
be obtained both in the direction of considering non-standard growth
conditions (see for instance the recent papers [23], [13], [15]), for which a
counterpart of [24] would be needed in the nonstandard setting, see [33]) and
in the direction of weakening the assumptions on the obstacle and on the
partial map $x\mapsto D_{\xi}f(x,\xi)$. We plan to continue our research in
these directions.

##  2\. Notations and statement of the main result

For a domain $\Omega\subset\mathbb{R}^{n}$, $n\geq 2$, and functions
$u:\Omega\rightarrow\mathbb{R},$ we consider the variational integral

$\mathcal{F}(u):=\int_{\Omega}f(x,Du)\,dx,$ (2.1)

where $f:\Omega\times\mathbb{R}^{n}\rightarrow[0,+\infty)$ is a Carathéodory
function which is convex and of class $\mathcal{C}^{2}$ with respect to the
second variable. Moreover, we suppose that there exist two positive constants
$\lambda,\Lambda$ such that, for all $\mu,\xi\in\mathbb{R}^{n}$,
$\mu=\mu_{i},$ $\xi=\xi_{i}$, $i=1,2,\dots,n$, a.e. in $\Omega$ and for $p\geq
2$

$\displaystyle\lambda\,(1+|\xi|^{2})^{\frac{p-2}{2}}|\mu|^{2}\leq\sum_{i,j}f_{\xi_{i}\xi_{j}}(x,\xi)\mu_{i}\mu_{j},$
(2.2)
$\displaystyle|f_{\xi_{i}\xi_{j}}(x,\xi)|\leq\,\Lambda\,(1+|\xi|^{2})^{\frac{p-2}{2}},$
(2.3) $\displaystyle\,\,|f_{\xi
x}(x,\xi)|\leq\,\Lambda\,(1+|\xi|^{2})^{\frac{p-1}{2}}.$ (2.4)

Throughout the paper we will denote by $B_{\rho}$ and $B_{R}$ balls of radii
respectively $\rho$ and $R$ (with $\rho<R$) compactly contained in $\Omega$
and with the same center, let us say $x_{0}\in\Omega$.

Our main result reads as follows

######  Theorem 2.1.

Let $u\in\mathcal{K}_{\psi,\Psi}(\Omega)$ be a solution to the obstacle
problem (1.1), under the assumptions (2.2)-(2.4). If $\psi\in
W^{2,\infty}_{\rm loc}(\Omega)$, then $u\in
W^{1,\infty}_{\textnormal{loc}}(\Omega)$ and the following estimate

$\|(1+|Du|^{2})^{\frac{1}{2}}\|_{L^{\infty}(B_{\rho})}\leq\,\frac{C}{(R-\rho)^{\beta}}\,\left(\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}}$
(2.5)

holds for every $0<\rho<R$ and with positive constants $C$ and $\beta$
depending on $n,p,\lambda,\Lambda,R,\rho$ and the local bounds for
$\|D\psi\|_{W^{1,\infty}}$.

######  Remark 2.2.

Due to the local nature of our main result, we do not assume further
regularity on the boundary datum $\Psi$.

##  3\. Proof of Theorem 2.1: the linearization procedure

Consider a smooth function $h_{\varepsilon}:(0,\infty)\to[0,1]$ such that
$h^{\prime}(s)\leq 0$ for all $s\in(0,\infty)$ and

$h(s)=\begin{cases}1&\quad\text{for}\quad s\leq\varepsilon,\\\
0&\quad\text{for}\quad s\geq 2\varepsilon.\\\ \end{cases}$

We remark that $u\in W^{1,p}(\Omega)$ is a solution to the obstacle problem in
$\mathcal{K}_{\psi,\Psi}(\Omega)$ if and only if $u\in\mathcal{K}_{\psi,\Psi}$
solves the following variational inequality

$\int_{\Omega}D_{\xi}f(x,Du)\cdot D(\varphi-u)\,dx\geq 0,$ (3.1)

for all $\varphi\in\mathcal{K}_{\psi,\Psi}$. In particular the function

$\varphi=u+t\cdot\eta\cdot h_{\varepsilon}(u-\psi)$

with $\eta\in C_{0}^{1}(\Omega)$, $\eta\geq 0$ and $0<t<<1$, is such that
$\varphi\in\mathcal{K}_{\psi,\Psi}$. Using $\varphi$ as test function, in the
variational inequality (3.1) gives

$\int_{\Omega}D_{\xi}f(x,Du)\cdot D(\eta h_{\varepsilon}(u-\psi))\,dx\geq
0\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$

Since

$\eta\mapsto L(\eta)=\int_{\Omega}D_{\xi}f(x,Du)\cdot D(\eta
h_{\varepsilon}(u-\psi))\,dx$

is a bounded positive linear functional, by the Riesz representation theorem
there exists a nonnegative measure $\lambda_{\varepsilon}$ such that

$\int_{\Omega}D_{\xi}f(x,Du)\cdot D(\eta
h_{\varepsilon}(u-\psi))\,dx=\int_{\Omega}\eta
d\lambda_{\varepsilon}\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$

We show that the measure $\lambda_{\varepsilon}$ is independent from
$\varepsilon$: indeed let $\varepsilon<\varepsilon^{\prime}$. The variation

$\varphi_{+}=u+t\cdot\eta\cdot(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi)$

is admissibile for all $\eta\in C_{0}^{1}(\Omega)$ and $0<t<<1$. This gives

$\int_{\Omega}D_{\xi}f(x,Du)\cdot
D(\eta(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi))\,dx\geq
0\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$

It’s easy to see that $(h_{\varepsilon}-h_{\varepsilon^{\prime}})(s)$ is
supported on $[\varepsilon,2\varepsilon^{\prime}]$, which implies

$\int_{\\{\varepsilon<u-\psi<2\varepsilon^{\prime}\\}}D_{\xi}f(x,Du)\cdot
D(\eta(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi))\,dx\geq
0\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$

Then, for $t<\frac{\varepsilon}{2},$ also the variation

$\varphi_{-}=u-t\cdot\eta\cdot(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi)$

is an admissibile test function; by the same argument we get

$\int_{\\{\varepsilon<u-\psi<2\varepsilon^{\prime}\\}}D_{\xi}f(x,Du)\cdot
D(\eta(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi))\,dx\leq
0\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,;$

together with the previous result, this implies

$\int_{\Omega}D_{\xi}f(x,Du)\cdot
D(\eta(h_{\varepsilon}-h_{\varepsilon^{\prime}})(u-\psi))\,dx=0\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,,$

that is equivalent to

$\int_{\Omega}\eta d\lambda_{\varepsilon}=\int_{\Omega}\eta
d\lambda_{\varepsilon^{\prime}}\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$
(3.2)

Now we can obtain equality between measures by a standard approximation
argument, let $\phi\in C_{0}^{\infty}(\Omega)$ be a smooth kernel with
$\phi\geq 0$ and $\int_{\Omega}\phi\,dx=1$, consider the corresponding family
of mollifiers $(\phi_{s})_{s>0}$ and any compact set $K$ such that
$K\subset\subset\Omega$. For every $s>0$ such that $s<\text{dist}(K,\Omega)$
we take in equation (3.2)

$\eta=\chi_{K}*\phi_{s}$

letting $s\searrow 0$ gives

$\int_{K}d\lambda_{\varepsilon}=\int_{K}d\lambda_{\varepsilon^{\prime}}\,,$

in conclusion $\lambda_{\varepsilon}=\lambda_{\varepsilon^{\prime}}$.
At this point we can rewrite our representation equation without the
$\varepsilon$ dependence on the measure

$\int_{\Omega}D_{\xi}f(x,Du)\cdot D(\eta
h_{\varepsilon}(u-\psi))\,dx=\int_{\Omega}\eta\,d\lambda\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,.$

By Theorem 2.1 in [24], that holds in particular under the assumptions
(2.2)-(2.4), we have that

$V_{p}(Du):=(1+|Du|^{2})^{\frac{p-2}{2}}Du\in W^{1,2}_{\rm loc}(\Omega),$

therefore we can integrate by parts and get

$\int_{\Omega}-\text{div}(D_{\xi}f(x,Du))\eta
h_{\varepsilon}(u-\psi)\,dx=\int_{\Omega}\eta\,d\lambda\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,.$

Now, in order to identify the measure $\lambda$, we may pass to the limit as
$\varepsilon\searrow 0$

$\int_{\Omega}-\text{div}(D_{\xi}f(x,Du))\chi_{\left[u=\psi\right]}\eta\,dx=\int_{\Omega}\eta\,d\lambda\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,.$ (3.3)

Let us introduce

$g:=\text{div}(D_{\xi}f(x,Du))\chi_{\left[u=\psi\right]};$ (3.4)

then the equation above takes the form:

$-\int_{\Omega}g\eta\,dx=\int_{\Omega}\eta\,d\lambda\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,.$

Now observing that

$\int_{\Omega}D_{\xi}f(x,Du)\cdot
D(\eta(1-h_{\varepsilon})(u-\psi))\,dx=0\qquad\forall\,\eta\in
C_{0}^{1}(\Omega)\,$

since $(1-h_{\varepsilon})(s)$ has support $[\varepsilon,+\infty)$, combining
our results we get

$\int_{\Omega}D_{\xi}f(x,Du)\cdot
D\eta\,dx=-\int_{\Omega}g\eta\,dx\qquad\forall\,\eta\in C_{0}^{1}(\Omega)\,.$
(3.5)

We are left to obtain an $L^{\infty}$ estimate for $g$: since $Du=D\psi$ a.e.
on the contact set, by (2.3) and (2.4) and the assumption $D\psi\in
W^{1,\infty}_{\rm loc}(\Omega;\mathbb{R}^{n})$, we have

$\begin{split}|g|&=\left|\text{div}(D_{\xi}f(x,Du))\chi_{\left[u=\psi\right]}\right|=\left|\text{div}(D_{\xi}f(x,D\psi))\right|\\\
&\leq\sum_{k=1}^{n}|f_{\xi_{k}x_{k}}(x,D\psi)|+\sum_{k,i=1}^{n}|f_{\xi_{k}\xi_{i}}(x,D\psi)\psi_{x_{k}x_{i}}|\\\
&\leq\Lambda(1+|D\psi|^{2})^{\frac{p-1}{2}}+\Lambda(1+|D\psi|^{2})^{\frac{p-2}{2}}|D^{2}\psi|\,\end{split}$

that is $g\in L^{\infty}_{\rm loc}(\Omega)$.

##  4\. Proof of Theorem 2.1: a priori estimate

Our starting point is now (3.5). We make use of the supplementary assumption
$u\in W^{1,\infty}_{\rm loc}(\Omega)$, which is needed in order to let (3.5)
to be satisfied; this assumption will be removed by means of the approximation
procedure. By the standard techniques of the difference quotients (see for
instance [35], [36]) we get

$u\in W^{2,2}_{\rm loc}(\Omega)\cap(1+|Du|^{2})^{\frac{p-2}{2}}|D^{2}u|^{2}\in
L^{1}_{\rm loc}(\Omega)$ (4.1)

so that the “second variation” system holds

$\int_{\Omega}\left(\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}D_{x_{i}}\varphi\,dx+\sum_{i=1}^{n}f_{\xi_{i}x_{s}}(x,Du)D_{x_{i}}\varphi\right)=\int_{\Omega}g\,D_{x_{s}}\varphi\,dx,$
(4.2)

for all $s=1,\dots,n$ and for all $\varphi\in W^{1,p}_{0}(\Omega)$; here $g$
is the function which has been introduced in (3.4). We fix $0<\rho<R$ with
$B_{R}$ compactly contained in $\Omega$ and we choose
$\eta\in\mathcal{C}^{1}_{0}(\Omega)$ such that $0\leq\eta\leq 1$, $\eta\equiv
1$ on $B_{\rho}$, $\eta\equiv 0$ outside $B_{R}$ and
$|D\eta|\leq\frac{C}{(R-\rho)}.$ We test (4.2) with
$\varphi=\eta^{2}(1+|Du|^{2})^{\gamma}u_{x_{s}}$, for some $\gamma\geq 0$
(which is possible due to our a priori assumption (4.1)) so that

$D_{x_{i}}\varphi=2\eta\eta_{x_{i}}(1+|Du|^{2})^{\gamma}u_{x_{s}}+2\eta^{2}\gamma(1+|Du|^{2})^{\gamma-1}|Du|D_{x_{i}}(|Du|)u_{x_{s}}+\eta^{2}(1+|Du|^{2})^{\gamma}u_{x_{s}x_{i}}$

Inserting in (4.2) we get:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}2\eta\eta_{x_{i}}(1+|Du|^{2})^{\gamma}u_{x_{s}}\,dx$
$\displaystyle+\int_{\Omega}\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}\eta^{2}(1+|Du|^{2})^{\gamma}u_{x_{s}x_{i}}\,dx$
$\displaystyle+\int_{\Omega}\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}2\eta^{2}\gamma(1+|Du|^{2})^{\gamma-1}|Du|D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{n}f_{\xi_{i}x_{s}}(x,Du)2\eta\eta_{x_{i}}(1+|Du|^{2})^{\gamma}u_{x_{s}}\,dx$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{n}f_{\xi_{i}x_{s}}(x,Du)\eta^{2}(1+|Du|^{2})^{\gamma}u_{x_{s}x_{i}}\,dx$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{n}f_{\xi_{i}x_{s}}(x,Du)2\eta^{2}\gamma(1+|Du|^{2})^{\gamma-1}|Du|D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle-\int_{\Omega}g2\eta\eta_{x_{s}}\,(1+|Du|^{2})^{\gamma}\,u_{x_{s}}\,dx$
$\displaystyle-\int_{\Omega}g2\eta^{2}\gamma\,(1+|Du|^{2})^{\gamma-1}|Du|\,D_{x_{s}}(|Du|)\,u_{x_{s}}\,dx$
$\displaystyle-\int_{\Omega}g\,\eta^{2}(1+|Du|^{2})^{\gamma}u_{x_{s}x_{s}}\,dx$
$\displaystyle=:I_{1,s}+I_{2,s}+I_{3,s}+I_{4,s}+I_{5,s}+I_{6,s}+I_{7,s}+I_{8,s}+I_{9,s}.$

We sum in the previous equation all terms with respect to $s$ from 1 to $n$,
and we denote by $I_{1}-I_{9}$ the corresponding integrals. In the sequel
constants will be denoted by $C$, regardless their actual value. Only the
relevant dependencies will be highlighted.
By the Cauchy-Schwarz inequality, the Young inequality and (2.3), we have

$\displaystyle|{I}_{1}|$ $\displaystyle=$
$\displaystyle\left|\int_{\Omega}2\eta(1+|Du|^{2})^{\gamma}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}\eta_{x_{i}}u_{x_{s}}\,dx\right|$
$\displaystyle=$
$\displaystyle\Bigg{|}\int_{\Omega}2\eta(1+|Du|^{2})^{\gamma}$
$\displaystyle\times\sum_{s=1}^{n}\left(\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)\eta_{x_{i}}\eta_{x_{j}}u_{x_{s}}^{2}\right)^{1/2}\left(\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{s}x_{i}}\,u_{x_{s}x_{j}}\right)^{1/2}\,dx\Bigg{|}$
$\displaystyle\leq$ $\displaystyle\int_{\Omega}2\eta(1+|Du|^{2})^{\gamma}$
$\displaystyle\times\left\\{\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)\eta_{x_{i}}\eta_{x_{j}}u_{x_{s}}^{2}\right\\}^{1/2}\,\left\\{\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{s}x_{i}}\,u_{x_{s}x_{j}}\right\\}^{1/2}\,dx$
$\displaystyle\leq$ $\displaystyle
C\int_{\Omega}(1+|Du|^{2})^{\gamma}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)\eta_{x_{i}}\eta_{x_{j}}u_{x_{s}}^{2}\,dx$
$\displaystyle+\frac{1}{4}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{s}x_{i}}\,u_{x_{s}x_{j}}\,dx$
$\displaystyle\leq$ $\displaystyle
C(\Lambda)\int_{\Omega}\Lambda(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}\sum_{i,j,s=1}^{n}\eta_{x_{i}}\eta_{x_{j}}u_{x_{s}}^{2}\,dx$
$\displaystyle+\frac{1}{4}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{s}x_{i}}\,u_{x_{s}x_{j}}\,dx$
$\displaystyle\leq$ $\displaystyle
C(n,\Lambda)\int_{\Omega}|D\eta|^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}\sum_{s=1}^{n}u_{x_{s}}^{2}\,dx$
$\displaystyle+\frac{1}{4}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{s}x_{i}}\,u_{x_{s}x_{j}}\,dx$
$\displaystyle\leq$ $\displaystyle
C\int_{\Omega}|D\eta|^{2}\,(1+|Du|^{2})^{p/2+\gamma}\,dx$
$\displaystyle+\frac{1}{4}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma}\,\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}u_{x_{i}x_{s}}\,dx.$

On the other hand, using (2.2) and the fact that
$D_{x_{j}}(|Du|)|Du|=\sum_{k=1}^{n}u_{x_{j}x_{k}}u_{x_{k}}$, we can estimate
the term $I_{3}$ as follows:

$\displaystyle|{I}_{3}|$ $\displaystyle=$
$\displaystyle\int_{\Omega}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)u_{x_{j}x_{s}}\left[2\eta^{2}\gamma(1+|Du|^{2})^{\gamma-1}D_{x_{i}}(|Du|)|Du|\right]u_{x_{s}}\,dx$
$\displaystyle\geq$ $\displaystyle
2\gamma\,\int_{\Omega}\eta^{2}|Du|^{2\gamma-1}\sum_{i,j,s=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)D_{x_{i}}(|Du|)u_{x_{j}x_{s}}u_{x_{s}}\,dx$
$\displaystyle=$ $\displaystyle
2\gamma\,\int_{\Omega}\eta^{2}|Du|^{2\gamma-1}\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)D_{x_{i}}(|Du|)\left(\sum_{s=1}^{n}u_{x_{j}x_{s}}u_{x_{s}}\right)\,dx$
$\displaystyle=$ $\displaystyle
2\gamma\,\int_{\Omega}\eta^{2}|Du|^{2\gamma}\sum_{i,j=1}^{n}f_{\xi_{i}\xi_{j}}(x,Du)D_{x_{i}}(|Du|)D_{x_{j}}(|Du|)\,dx$
$\displaystyle\geq$ $\displaystyle
2\gamma\,\int_{\Omega}\eta^{2}|Du|^{2\gamma}|D(|Du|)|^{2}(1+|Du|^{2})^{\frac{p-2}{2}}\,dx\geq
0.$

We can estimate the fourth and the fifth term by the Cauchy-Schwarz and the
Young inequalities, together with (2.4), as follows

$\displaystyle|{I}_{4}|$ $\displaystyle=$ $\displaystyle
2\int_{\Omega}\eta(1+|Du|^{2})^{\gamma}\sum_{i,s=1}^{n}f_{\xi_{i}x_{s}}(x,Du)\eta_{x_{i}}u_{x_{s}}\,dx$
$\displaystyle\overset{\eqref{gtx}}{\leq}$ $\displaystyle
2\Lambda\,\int_{\Omega}\eta\,(1+|Du|^{2})^{\gamma+\frac{p-1}{2}}\sum_{i,s=1}^{n}|\eta_{x_{i}}u_{x_{s}}|\,dx$
$\displaystyle\leq$ $\displaystyle
C(\Lambda)\int_{\Omega}\eta|D\eta||Du|\,(1+|Du|^{2})^{\gamma+\frac{p-1}{2}}\,dx$
$\displaystyle\leq$ $\displaystyle
C\int_{\Omega}(\eta^{2}+|D\eta|^{2})(1+|Du|^{2})^{\gamma+\frac{p}{2}},$

and also:

$\displaystyle|{I}_{5}|$ $\displaystyle=$
$\displaystyle\left|\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma}\sum_{i,s=1}^{n}f_{\xi_{i}x_{s}}(x,Du)u_{x_{s}x_{i}}\,dx\right|$
$\displaystyle\overset{\eqref{gtx}}{\leq}$
$\displaystyle\Lambda\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma+\frac{p-1}{2}}\,\sum_{i,s=1}^{n}u_{x_{s}x_{i}}\,dx$
$\displaystyle=$
$\displaystyle\Lambda\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma+\frac{p-1}{2}}|D^{2}u|\,dx$
$\displaystyle\leq$
$\displaystyle\frac{1}{4}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}|D^{2}u|^{2}\,dx+C\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx.$

Finally, by similar arguments, the Cauchy-Schwarz inequality, (2.3) and
$|D(|Du|)|\leq|D^{2}u|$ give

$\displaystyle|{I}_{6}|$ $\displaystyle=$ $\displaystyle
2\gamma\,\int_{\Omega}\sum_{i,s=1}^{n}f_{\xi_{i}x_{s}}(x,Du)\eta^{2}(1+|Du|^{2})^{\gamma-1}|Du|D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle=$ $\displaystyle
2\gamma\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma-1}|Du|\sum_{i,s=1}^{n}f_{\xi_{i}x_{s}}(x,Du)D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle\leq$ $\displaystyle
2\,\gamma\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma-\frac{1}{2}}\sum_{i,s=1}^{n}f_{\xi_{i}x_{s}}(x,Du)D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle\leq$ $\displaystyle
2\,\gamma\,\Lambda\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma-\frac{1}{2}+\frac{p-1}{2}}\sum_{i,s=1}^{n}D_{x_{i}}(|Du|)u_{x_{s}}\,dx$
$\displaystyle\leq$ $\displaystyle
2\,n\,\gamma\,\Lambda\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\gamma-\frac{1}{2}+\frac{p-1}{2}}|D(|Du|)||Du|\,dx$
$\displaystyle\leq$ $\displaystyle
2\,n\,\gamma\,\Lambda\,\int_{\Omega}\,\eta^{2}(1+|Du|^{2})^{\gamma+\frac{p-1}{2}}|D^{2}u|\,dx$
$\displaystyle\leq$
$\displaystyle\frac{1}{4}\int_{\Omega}\eta^{2}|D^{2}u|^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}\,dx+C\gamma^{2}\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx,$

where the constant $C$ depends only on $n,p,\Lambda$ but it is independent of
$\gamma$.
Let us now deal with the terms containing the function $g$. We use the bound

$\|g\|_{L^{\infty}_{\rm loc}(\Omega)}\leq\,C,$

established in Section 3; here the constant $C$ is independent of $n$ and
$\gamma$.
We first have

$\displaystyle|I_{7}|$ $\displaystyle\leq$
$\displaystyle\|g\|_{L^{\infty}(B_{R})}\int_{\Omega}2|D\eta|(1+|Du|^{2})^{\gamma}\,|Du|\,dx$
$\displaystyle\leq$
$\displaystyle\,C\,\int_{\Omega}|D\eta|^{2}(1+|Du|^{2})^{\gamma+p/2}\,dx.$

Arguing similarly as in the estimate of $I_{6}$ we deduce

$\displaystyle|I_{8}|$ $\displaystyle\leq$
$\displaystyle\,2\gamma\int_{\Omega}|g|\,\eta^{2}(1+|Du|^{2})^{\gamma-1}|D(|Du|)|\,|Du|^{2}\,dx$
$\displaystyle\leq$
$\displaystyle\frac{1}{4}\int_{\Omega}\eta^{2}|D^{2}u|^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}\,dx+C\gamma^{2}\,\|g\|^{2}_{L^{\infty}(B_{R})}\,\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx.$

Finally

$\displaystyle|I_{9}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{4}\int_{\Omega}\eta^{2}|D^{2}u|^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}\,dx+C\|g\|^{2}_{L^{\infty}(B_{R})}\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx.$

We remark that in the estimates of the terms $I_{7}-I_{9}$ we used in an
essential way the fact $(1+|Du|^{2})\geq 1$. It is possible to overcome this
limtation and dealing also with the degenerate case by proceeding as in [23].
Summing up and using (2.2) we obtain

$\int_{\Omega}\,\eta^{2}\,(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}|D^{2}u|^{2}\,dx\leq\,C\,(1+\gamma^{2})\int_{\Omega}(\eta^{2}+|D\eta|^{2})\,(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx,$

for any $0<\rho<R$, where the constant $C$ depends on $\lambda,\Lambda,n,p$
but is independent of $\gamma$.

On the other hand

$\displaystyle\int_{\Omega}\left|D\left[\eta(1+|Du|^{2})^{\frac{p}{4}+\frac{\gamma}{2}}\right]\right|^{2}\,dx$
$\displaystyle\leq$
$\displaystyle\int_{\Omega}|D\eta|^{2}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx+\int_{\Omega}\eta^{2}(1+|Du|^{2})^{\frac{p-2}{2}+\gamma}|D^{2}u|\,dx$
$\displaystyle\leq$
$\displaystyle\,C\,(1+\gamma^{2})\int_{\Omega}(\eta^{2}+|D\eta|^{2})\,(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx.$

Now, let $2^{*}=\frac{2n}{n-2}$ for $n>2$ while $2^{*}$ equal to any fixed
real number greater than 2, if $n=2$. By the Sobolev’s inequality there exists
a constant $C$ such that

$\displaystyle\left(\int_{\Omega}\eta^{2^{*}}\left(1+|Du|^{2}\right)^{\left(\frac{p}{2}+\gamma\right)\frac{2^{*}}{2}}\,dx\right)^{\frac{2}{2^{*}}}$
$\displaystyle\leq$
$\displaystyle\,C\,\int_{\Omega}\left|D\left[\eta(1+|Du|^{2})^{\frac{p}{4}+\frac{\gamma}{2}}\right]\right|^{2}\,dx$
$\displaystyle\leq$
$\displaystyle\,C\,(1+\gamma^{2})\int_{\Omega}(\eta^{2}+|D\eta|^{2})\,(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx.$

##  5\. Proof of Theorem 2.1: an iteration scheme

By the definition of $\eta$ we obtain

$\left(\int_{B_{\rho}}\left(1+|Du|^{2}\right)^{\left(\frac{p}{2}+\gamma\right)\frac{2^{*}}{2}}\,dx\right)^{\frac{2}{2^{*}}}\leq\frac{C(1+\gamma^{2})}{(R-\rho)^{2}}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}+\gamma}\,dx,$
(5.1)

for any $B_{\rho}$ and $B_{R}$ compactly contained in $\Omega$ with $\rho<R$
and the same center.
Let us define the sequence $(\gamma_{k})_{k}$ by induction as follows

$\gamma_{1}=0\,;\qquad\gamma_{k+1}=\left(\gamma_{k}+\frac{p}{2}\right)\frac{2^{*}}{2}-\frac{p}{2},\quad\forall\,k\geq
1.$ (5.2)

If $(\gamma_{k})_{k}$ is the sequence defined in (5.2), then the following
representation formula holds

$\gamma_{k}=\frac{p}{2}\left(\left(\frac{2^{*}}{2}\right)^{k-1}-1\right),$
(5.3)

the proof easily follows by induction.

Fix two rays $R_{0}$ and $\rho_{0}$ such that $0<\rho_{0}<R_{0}<\rho_{0}+1$,
let us define $R_{k}=\rho_{0}+(R_{0}-\rho_{0})2^{-k}$ for every $k\geq 1$.
Since $R_{k+1}<R_{k}\,\,\forall\,k\geq 1$, we are allowed to insert
$\rho=R_{k+1}$ and $R=R_{k}$ in estimate (5.1), we also take
$\gamma=\gamma_{k}$, where $(\gamma_{k})_{k}$ has been defined in (5.2). With
these choices, by (5.1) we obtain

$\displaystyle\left(\int_{B_{R_{k+1}}}(1+|Du|^{2})^{\left(\frac{p}{2}+\gamma_{k}\right)\frac{2^{*}}{2}}\,dx\right)^{\frac{2}{2^{*}}}$
$\displaystyle\leq\frac{C(1+\gamma_{k}^{2})}{(R_{k}-R_{k+1})^{2}}\int_{B_{R_{k}}}(1+|Du|^{2})^{\frac{p}{2}+\gamma_{k}}\,dx$
$\displaystyle=\frac{C(1+\gamma_{k}^{2})4^{k+1}}{(R_{0}-\rho_{0})^{2}}\int_{B_{R_{k}}}(1+|Du|^{2})^{\frac{p}{2}+\gamma_{k}}\,dx,$
(5.4)

for every $k\geq 1$. Since the constant $C$ was independent of $\gamma$, it is
clear that now $C$ is independent on $k$. Let us define

$A_{k}=\left(\int_{B_{R_{k}}}(1+|Du|^{2})^{\frac{p}{2}+\gamma_{k}}\,dx\right)^{\frac{1}{p+2\gamma_{k}}},\quad\forall\,k\geq
1.$ (5.5)

The estimate (5.4) directly gives us an inductive estimate on $(A_{k})_{k}$.
By definition (5.2) of $\gamma_{k}$, and by (5.4) we have

$\displaystyle A_{k+1}$
$\displaystyle=\left(\int_{B_{R_{k+1}}}(1+|Du|^{2})^{\frac{p}{2}+\gamma_{k+1}}\,dx\right)^{\frac{1}{p+2\gamma_{k+1}}}$
$\displaystyle=\left(\int_{B_{R_{k+1}}}(1+|Du|^{2})^{\left(\frac{p}{2}+\gamma_{k}\right)\frac{2^{*}}{2}}\,dx\right)^{\frac{1}{\left(\frac{p}{2}+\gamma_{k}\right)2^{*}}}$
$\displaystyle\leq\left[\frac{C(1+\gamma_{k}^{2})4^{k+1}}{(R_{0}-\rho_{0})^{2}}\right]^{\frac{1}{p+2\gamma_{k}}}\cdot\left(\int_{B_{R_{k}}}(1+|Du|^{2})^{\frac{p}{2}+\gamma_{k}}\,dx\right)^{\frac{1}{p+2\gamma_{k}}}$
$\displaystyle=\left[\frac{C(1+\gamma_{k}^{2})4^{k+1}}{(R_{0}-\rho_{0})^{2}}\right]^{\frac{1}{p+2\gamma_{k}}}\cdot
A_{k}\,,$ (5.6)

for every $k\geq 1$. By iterating (5.6) we obtain

$A_{k+1}\leq
A_{1}\cdot\prod_{i=1}^{k}\left[\frac{C(1+\gamma_{i}^{2})4^{i+1}}{(R_{0}-\rho_{0})^{2}}\right]^{\frac{1}{p+2\gamma_{i}}}\leq
C\cdot A_{1}\cdot(R_{0}-\rho_{0})^{-\sum_{i=1}^{\infty}1/(p/2+\gamma_{i})},$
(5.7)

for every $k\geq 1$. Here we used the fact that
$\prod_{i=1}^{k}\big{[}C(1+\gamma_{i}^{2})4^{i+1}\big{]}^{\frac{1}{p+2\gamma_{i}}}$
is bounded, indeed

$\prod_{i=1}^{\infty}\big{[}C(1+\gamma_{i}^{2})4^{i+1}\big{]}^{\frac{1}{p+2\gamma_{i}}}=\textnormal{exp}\left(\sum_{i=1}^{\infty}\frac{\ln(C(1+\gamma_{i}^{2})4^{i+1})}{p+2\gamma_{i}}\right)<+\infty\,,$

since the series is convergent by formula (5.3). Since $\rho_{0}<R_{k}<R_{0}$,
by exploting $A_{1}$ in (5.7), for every $k\geq 2$ we deduce that

$\left(\int_{B_{\rho_{0}}}(1+|Du|^{2})^{\frac{1}{2}\left(p+2\gamma_{k}\right)}\,dx\right)^{\frac{1}{p+2\gamma_{k}}}\leq\frac{C}{(R_{0}-\rho_{0})^{\beta}}\left(\int_{B_{R_{0}}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}},$
(5.8)

where we defined

$\beta:=\sum_{i=1}^{\infty}\frac{1}{\frac{p}{2}+\gamma_{i}}=\frac{2^{*}}{p(2^{*}/2)-p}.$

Since $p+2\gamma_{k}\to+\infty$, as $k\to+\infty$ it is clear that the left
hand side of (5.8) converges to the essential supremum of
$(1+|Du|^{2})^{\frac{1}{2}}$ in $B_{\rho_{0}}$. Thus, provided the further
regulatity that we will remove in the next section with an approximation
procedure, we conclude the proof of Theorem 2.1.

##  6\. Proof of Theorem 2.1: approximation

First of all we state an approximation theorem for $f$ through a suitable
sequence of regular functions. Let $B$ be the unit ball of $\mathbb{R}^{n}$
centered in the origin and consider a positive decreasing sequence
$\varepsilon_{\ell}\rightarrow 0$. Let $\phi\in C_{0}^{\infty}(B)$ be a smooth
symmetric kernel with $\phi\geq 0$ and $\int_{B}\phi\,dx=1$.
Let us define

$f^{\ell}(x,\xi)=\int_{B\times
B}\phi(y)\phi(\eta)f(x+\varepsilon_{\ell}y,\xi+\varepsilon_{\ell}\eta)\,d\eta\,dy,$
(6.1)

and consequently

$f^{\ell k}(x,\xi)=f^{\ell}(x,\xi)+\frac{1}{k}(1+|\xi|^{2})^{\frac{p}{2}}.$
(6.2)

The proof of the following theorem follows from similar arguments as in [23].

######  Theorem 6.1.

Let $f$ satisfy the growth conditions (2.2), (2.3), (2.4), and let $f(x,\xi)$
be a convex $C^{2}$ function in its second argument. Then the sequence
$f^{\ell k}:\Omega\times\mathbb{R}^{n}\rightarrow[0,+\infty)$ defined in
(6.2), convex in the last variable is such that $f^{\ell k}\to f$ pointwise as
$(\ell,k)\to+\infty$ for a.e. $x\in\Omega$ and for all $\xi\in\mathbb{R}^{n}$.
Moreover, $f^{\ell k}\to f$ uniformly in $\Omega_{0}\times K$ for every
$\Omega_{0}\subset\subset\Omega$ and $K$ compact set of $\mathbb{R}^{n}$.
Moreover:

  * •

There exist constants $M_{1},M_{2}>0$, independent of $k,\ell$ such that

$M_{1}(1+|\xi|^{2})^{\frac{p}{2}}\leq f^{\ell k}(x,\xi)\leq
M_{2}(1+|\xi|^{2})^{\frac{p}{2}}\qquad\textnormal{for a.e. $x\in\Omega$, for
all $\xi\in\mathbb{R}^{n}$},$ (6.3)

  * •

there exist constants $\lambda_{1},\Lambda_{1}>0$ such that for all
$(x,\xi)\in\Omega\times\mathbb{R}^{n}$ and $\mu\in\mathbb{R}^{n}$

$\lambda_{1}\,(1+|\xi|^{2})^{\frac{p-2}{2}}|\mu|^{2}\leq\sum_{i,j}f^{\ell
k}_{\xi_{i}\xi_{j}}(x,\xi)\mu_{i}\mu_{j},$ (6.4)

and

$|f^{\ell
k}_{\xi_{i}\xi_{j}}(x,\xi)|\leq\Lambda_{1}\,(1+|\xi|^{2})^{\frac{p-2}{2}}.$
(6.5)

  * •

There exists a constant $\Lambda_{2}>0$ such that for all
$(x,\xi)\in\Omega\times\mathbb{R}^{n}$

$|f^{\ell k}_{\xi x}(x,\xi)|\leq\Lambda_{2}(1+|\xi|^{2})^{\frac{p-1}{2}}.$
(6.6)

Let $u\in W^{1,p}(\Omega)$ denote a local minimizer of the problem

$\min\left\\{\int_{\Omega}f(x,Dw):w\in\mathcal{K}_{\psi,\Psi}(\Omega)\right\\}.$
(6.7)

Fix $x_{0}\in\Omega$ and a radius $R>0$ such that
$B_{R}(x_{0})\subset\subset\Omega$. Consider the following variational problem

$\inf\left\\{\int_{B_{R}}f^{\ell k}(x,Dv)\,dx,\,\,\,v\in
W^{1,p}_{0}(B_{R})+u\,,\,v\geq\psi\,\,\,\textnormal{a.e. in
$\Omega$}\,\right\\},$ (6.8)

where $f^{\ell k}$ have been defined in (6.2). By classical lower
semicontinuity arguments, there exist $v^{\ell k}\in W^{1,p}_{0}(B_{R})+u$
solution to problem (6.8). By the growth conditions (6.3) and the minimality
of $v^{\ell k}$, we get

$\displaystyle\int_{B_{R}}|Dv^{\ell k}|^{p}\,dx$
$\displaystyle\leq\int_{B_{R}}(1+|Dv^{\ell k}|^{2})^{\frac{p}{2}}\,dx$
$\displaystyle\leq C\int_{B_{R}}f^{\ell k}(x,Dv^{\ell k})\,dx\leq
C\int_{B_{R}}f^{\ell k}(x,Du)\,dx$
$\displaystyle=C\int_{B_{R}}f^{\ell}(x,Du)\,dx+\frac{C}{k}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx$
$\displaystyle\leq C\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx.$

Note that the right hand side of the previous estimate does not depend on
$\ell$, thus $\|Dv^{\ell k}\|_{L^{p}(B_{R})}$ is bounded for every fixed $k$.
Since $W^{1,p}(B_{R})$ is weakly compact, up to passing to a subsequence, we
can assume there exist $v^{k}\in u+W^{1,p}_{0}(B_{R})$ such that

$v^{\ell
k}\stackrel{{\scriptstyle\ell\rightarrow\infty}}{{\rightarrow}}v^{k}\,\,\,\textnormal{weakly
in $W^{1,p}_{0}(B_{R})+u$},$ (6.9)

where, for every $k\geq 1$, the limit function $v^{k}$ still belongs to
$\mathcal{K}_{\psi,\Psi}$ since this set is weakly closed. Moreover, by
classical convolution properties $f^{\ell}(x,Du)\to f(x,Du)$ pointwise a.e. in
$B_{R}$ as $\ell\to+\infty$, and again by growth conditions (6.3)

$\int_{B_{R}}f^{\ell}(x,Du)\,dx\leq
C\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx,$ (6.10)

for every $\ell\geq 1$. Thus by the Dominated Convergence Theorem we deduce

$\lim_{\ell\rightarrow\infty}\int_{B_{R}}|Dv^{\ell k}|^{p}\,dx\leq
C\int_{B_{R}}f(x,Du)\,dx+\frac{C}{k}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx.$
(6.11)

Now, in force of Theorem 6.1, we know that $f^{\ell k}$ satisfies the growth
conditions needed to apply our a priori estimate for every $\ell,k$. Thus,
applying the a priori estimate to the solutions to (6.8), using the minimality
of $v^{\ell k}$

$\displaystyle\|Dv^{\ell k}\|_{L^{\infty}(B_{\rho})}$
$\displaystyle\leq\widetilde{C}\left(\int_{B_{R}}(1+|Dv^{\ell
k}|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}}$
$\displaystyle\leq\widetilde{C}\left(\int_{B_{R}}f^{\ell k}(x,Dv^{\ell
k})\,dx\right)^{\frac{1}{p}}$
$\displaystyle\leq\widetilde{C}\left(\int_{B_{R}}f^{\ell
k}(x,Du)\,dx\right)^{\frac{1}{p}}$
$\displaystyle\leq\widetilde{C}\left(\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}},$

where the constant $\widetilde{C}=\frac{C}{(R-\rho)^{\beta}}$ is independent
of $\ell,k$. Therefore, for every fixed $k$, the sequence $\|Dv^{\ell
k}\|_{L^{\infty}(B_{R})}$ is bounded. Since weak limits are unique, we
conclude that

$v^{\ell
k}\stackrel{{\scriptstyle\ell\rightarrow\infty}}{{\rightarrow}}v^{k}\,\,\,\textnormal{weakly
star in $W^{1,\infty}_{\rm loc}(B_{R})$}.$ (6.12)

By estimate (6.11)

$\displaystyle\|Dv^{k}\|_{L^{p}(B_{R})}$
$\displaystyle\leq\liminf_{\ell\to\infty}\|Dv^{\ell k}\|_{L^{p}(B_{R})}$
$\displaystyle\leq
C\left(\int_{B_{R}}f(x,Du)\,dx+\frac{1}{k}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}}$
$\displaystyle\leq
C\left(\int_{B_{R}}f(x,Du)+(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}},$
(6.13)

for every $k\geq 1$. Moreover, for any $\rho$ such that $0<\rho<R$, we have

$\displaystyle\|Dv^{k}\|_{L^{\infty}(B_{\rho})}$
$\displaystyle\leq\liminf_{\ell\to\infty}\|Dv^{\ell
k}\|_{L^{\infty}(B_{\rho})}$
$\displaystyle\leq\widetilde{C}\left(\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}}.$
(6.14)

Thus we can deduce that there exists $\bar{v}\in u+W^{1,p}_{0}(B_{R})$ such
that, up to subsequences

$v^{k}\rightarrow\bar{v}\,\,\,\textnormal{weakly in
$W^{1,p}_{0}(B_{R})+u$},\qquad\qquad
v^{k}\rightarrow\bar{v}\,\,\,\textnormal{weakly star in $W^{1,\infty}_{\rm
loc}(B_{R})$}$

and the limit function $\bar{v}$ still belongs to $\mathcal{K}_{\psi,\Psi}$
since this set is weakly closed. Since estimate (6.14) holds for every $k\geq
1$, for any $0<\rho<R$ we obtain

$\|D\bar{v}\|_{L^{\infty}(B_{\rho})}\leq\frac{C}{(R-\rho)^{\beta}}\left(\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)^{\frac{1}{p}}.$
(6.15)

Hence combining lower semicontinuity, the minimality of $v^{\ell k}$ for
problem (6.8), and properties of mollification we obtain

$\displaystyle\int_{B_{\rho}}f(x,D\bar{v})\,dx$
$\displaystyle\leq\liminf_{k\to\infty}\int_{B_{\rho}}f(x,Dv^{k})\,dx$
$\displaystyle\leq\liminf_{k\to\infty}\liminf_{\ell\to\infty}\int_{B_{\rho}}f(x,Dv^{\ell
k})\,dx$
$\displaystyle=\liminf_{k\to\infty}\liminf_{\ell\to\infty}\int_{B_{\rho}}f^{\ell}(x,Dv^{\ell
k})\,dx$
$\displaystyle\leq\liminf_{k\to\infty}\liminf_{\ell\to\infty}\left(\int_{B_{R}}f^{\ell}(x,Dv^{\ell
k})+\frac{1}{k}(1+|Dv^{\ell k}|^{2})^{\frac{p}{2}}\,dx\right)$
$\displaystyle=\liminf_{k\to\infty}\liminf_{\ell\to\infty}\int_{B_{R}}f^{\ell
k}(x,Dv^{\ell k})\,dx$
$\displaystyle\leq\liminf_{k\to\infty}\liminf_{\ell\to\infty}\int_{B_{R}}f^{\ell
k}(x,Du)\,dx$
$\displaystyle=\liminf_{k\to\infty}\liminf_{\ell\to\infty}\left(\int_{B_{R}}f^{\ell}(x,Du)\,dx+\frac{1}{k}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)$
$\displaystyle=\liminf_{k\to\infty}\left(\int_{B_{R}}f(x,Du)\,dx+\frac{1}{k}\int_{B_{R}}(1+|Du|^{2})^{\frac{p}{2}}\,dx\right)$
$\displaystyle=\int_{B_{R}}f(x,Du)\,dx,$

for every $0<\rho<R$. Then, letting $\rho\rightarrow R$ in the previous
inequality, we get

$\int_{B_{R}}f(x,D\bar{v})\,dx\leq\int_{B_{R}}f(x,Du)\,dx.$ (6.16)

Therefore, $u$ and $\bar{v}$ are two solutions to problem (6.7). We claim that
$u=\bar{v}$ on $B_{R}$. Indeed, if we suppose that $u\neq{v}$ choose any
$\theta\in(0,1)$ and define $v_{\theta}=\theta\bar{v}+(1-\theta)u$, since
$f(x,\xi)$ is strictly convex in its second argument we have

$\displaystyle\int_{B_{R}}f(x,Dv_{\theta})\,dx$
$\displaystyle<\theta\int_{B_{R}}f(x,D\bar{v})\,dx+(1-\theta)\int_{B_{R}}f(x,Du)\,dx$
$\displaystyle=\int_{B_{R}}f(x,Du)\,dx,$

contradicting the minimality of $u$. Thus, regularity is preserved through the
limit procedure, and this concludes the proof of Theorem 2.1 because we have
removed the further regularity assumptions we needed to obtain the result.

## References

  * [1] C. Baiocchi, A. Capelo: _Disequazioni variazionali e quasi variazionali. Applicazioni a problemi di frontiera libera,_ Quaderni U.M.I. Pitagora, Bologna 1978. English transl. J. Wiley, Chichester-New York, 1984.
  * [2] P. Baroni: _Lorentz estimates for obstacle parabolic problems_ , Nonlinear Anal., 96, (2014), 167-188.
  * [3] M. Bildhauer, M. Fuchs, G. Mingione: _A priori gradient bounds and local $\mathcal{C}^{1,\alpha}$-estimates for (double) obstacle problems under non-standard growth conditions_, Z. Anal. Anwendungen, 20, (4), (2001), 959-985.
  * [4] V. Bögelein, F. Duzaar, G. Mingione: _Degenerate problems with irregular obstacles_ , J. Reine Angew. Math., 650, (2011), 107-160.
  * [5] V. Bögelein, T. Lukkari, C. Scheven: _The obstacle problem for the porous medium equation_ , Math. Ann. 363, (1) (2015), 455-499.
  * [6] H. Brézis, D. Kinderlehrer: _The smoothness of solutions to nonlinear variational inequalities_ , Indiana Univ. Math. J., 23 (1973-1974), 831-844.
  * [7] S.S. Byun, Y. Cho, J. Ok: _Global gradient estimates for nonlinear obstacle problem with non standard growth,_ Forum Math., 28, (4), (2016), 729-747.
  * [8] L. Caffarelli: _The regularity of elliptic and parabolic free boundaries_ , Bull. Amer. Math. Soc., 82, (1976), 616-618.
  * [9] L.A. Caffarelli, D. Kinderlehrer: _Potential methods in variational inequalities_ , J. Anal. Math., 37 (1980), 285-295.
  * [10] M. Chipot: _Variational inequalities and flow in porous media_ , Springer-Verlag, Berlin and New York, 1984.
  * [11] H. Choe: _A regularity theory for a general class of quasilinear elliptic partial differential equations and obstacle problems_ , Arch. Rational Mech. Anal., 114, (4), (1991), 383-394.
  * [12] H. Choe, J.-L. Lewis: _On the obstacle problem for quasilinear elliptic equations of $p$ Laplacian type_, SIAM J. Math. Anal., 22, (3), (1991), 623-638.
  * [13] G. Cupini, F. Giannetti, R. Giova, A. Passarelli di Napoli: _Higher integrability for minimizers of asymptotically convex integrals with discontinuous coefficients_ , Nonlinear Anal., 154, (2017), 7-24.
  * [14] G. Cupini, F. Giannetti, R. Giova, A. Passarelli di Napoli: _Regularity results for vectorial minimizers of a class of degenerate convex integrals,_ J. Differential Equations 265, (2018), no. 9, 4375–4416.
  * [15] G. Cupini, P. Marcellini, E. Mascolo: _Existence and regularity for elliptic equations under p,q -growth,_ Advances in Differential Equations, 19(7/8), (2014), 693-724.
  * [16] F. Duzaar: _Variational inequalities for harmonic maps,_ J. Reine Angew. Math., 374 (1987), 39-60.
  * [17] M. Eleuteri: Regularity results for a class of obstacle problems, Appl. Math., 52, (2), (2007), 137-169.
  * [18] M. Eleuteri, J. Habermann: Regularity results for a class of obstacle problems under non standard growth conditions, J. Math. Anal. Appl., 344, (2), (2008), 1120-1142.
  * [19] M. Eleuteri, J. Habermann: Calderón-Zygmund type estimates for a class of obstacle problems with $p(x)$ growth, J. Math. Anal. Appl., 372, (1), (2010), 140-161.
  * [20] M. Eleuteri, J. Habermann: A Hölder continuity result for a class of obstacle problems under non standard growth conditions, Math. Nachr., 284, No. 11-12, (2011), 1404-1434.
  * [21] M. Eleuteri, P. Harjulehto, T. Lukkari: Global regularity and stability of solutions to obstacle problems with nonstandard growth, Rev. Mat. Complut., 26, No. 1, (2013), 147-181.
  * [22] M. Eleuteri, P. Marcellini, E. Mascolo: Lipschitz estimates for systems with ellipticity conditions at infinity, Ann. Mat. Pura Appl., 195 (5), (2016), 1575-1603.
  * [23] M. Eleuteri, P. Marcellini, E. Mascolo: _Regularity for scalar integrals without structure conditions,_ Advances in Calculus of Variations, (2018), DOI: https://doi.org/10.1515/acv-2017-0037
  * [24] M. Eleuteri, A. Passarelli di Napoli: Higher differentiability for solutions to a class of obstacle problems, Calc. Var. Partial Differential Equations, 57: 115 (2018).
  * [25] M. Eleuteri, A. Passarelli di Napoli: Regularity results for a class of non-differentiable obstacle problems, Nonlinear Anal, (2019) (to appear). DOI: 10.1016/j.na.2019.01.024
  * [26] G. Fichera: _Problemi elastostatici con vincoli unilaterali: il problema di Signorini con ambigue condizioni al contorno,_ Atti Accad. Naz. Lincei Mem. Cl. Sci. Fis. Mat. Nat. Sez. Ia, 7, (8), (1963-1964), 91-140.
  * [27] A. Friedman: _Variational principles and free boundary problems_ , Wiley, New York, 1982.
  * [28] M. Fuchs: _Variational inequalities for vector valued functions with non convex obstacles,_ Analysis 5 (1985), 223-238.
  * [29] M. Fuchs: _p-harmonic obstacle problems. Part I: Partial regularity theory,_ Ann. Mat. Pura Appl. 156 (1990), 127-158.
  * [30] M. Fuchs: _Topics in the Calculus of Variations_ , Advanced Lectures in Mathematics, Vieweg, (1994).
  * [31] M. Fuchs, L. Gongbao: _Variational inequalities for energy functionals with nonstandard growth conditions_ , Abstr. Appl. Anal., 3 (1998), 41-64.
  * [32] M. Fuchs, G. Mingione: _Full $\mathcal{C}^{1,\alpha}-$regularity for free and constrained local minimizers of elliptic variational integrals with nearly linear growth_, Manuscripta Math., 102, (2000), 227-250.
  * [33] C. Gavioli: _Higher differentiability for a class of obstacle problems with nonstandard growth conditions_ , (2019) Forum Mathematicum. 2019, Vol. 31, No. 6, pp. 1501-1516.
  * [34] C. Gavioli: _A priori estimates for solutions to a class of obstacle problems under p, q- growth conditions._ Journal of Elliptic and Parabolic Equations. 2019, Vol. 5, No. 2, pp. 325-347.
  * [35] M. Giaquinta:_Multiple integrals in the calculus of variations and nonlinear elliptic systems,_ Annals of Mathematics Studies, 105. Princeton University Press. (1983).
  * [36] E. Giusti:_Direct methods in the calculus of variations,_ World Scientific Publishing Co., Inc., River Edge, NJ (2003).
  * [37] D. Kinderlehrer, G. Stampacchia: _An introduction to variational inequalities and their applications_ , Academic Press, NY (1980).
  * [38] P. Lindqvist: _Regularity for the gradient of the solution to a nonlinear obstacle problem with degenerate ellipticity_ , Nonlinear Anal., 12, (11), (1988), 1245-1255.
  * [39] J.L. Lions, G. Stampacchia: _Variational inequalities_ , Comm. Pure Appl. Math., 20, (1967), 493-519.
  * [40] J.H. Michael, W.P. Ziemer: _Interior regularity for solutions to obstacle problems_ , Nonlinear Anal., 10, (12), (1986), 1427-1448.
  * [41] P. Marcellini: Regularity of minimizers of integrals in the calculus of variations with non standard growth conditions,_Arch. Rational Mech. Anal._ 105 (1989) 267-284.
  * [42] P. Marcellini: Regularity and existence of solutions of elliptic equations with $p-q$-growth conditions,_J. Differential Equations_ 90 (1991) 1-30.
  * [43] P. Marcellini: Regularity for elliptic equations with general growth conditions,_J. Differential Equations_ 105 (1993) 296-333.
  * [44] P. Marcellini: Regularity for some scalar variational problems under general growth conditions,_J. Optim. Theory Appl._ 90, no. 1 (1996) 161-181.
  * [45] J. Moser: _A new proof of De Giorgi’s theorem concerning the regularity problem for elliptic differential equations,_ Comm. Pure Appl. Math., 13 (1960), 475-468.
  * [46] J. Ok: _Regularity results for a class of obstacle problems with nonstandard growth_ , J. Math. Anal. Appl., 444, (2), (2016), 957-979.
  * [47] J. Ok: _Calderón-Zygmund estimates for a class of obstacle problems with nonstandard growth_ , NoDEA, Nonlinear Diff. Equ. Appl., 23, (4), (2016).
  * [48] J. Ok: _Gradient continuity for nonlinear obstacle problems_ , Mediterr. J. Math., 14, (1), (2017).
  * [49] J.F. Rodrigues: _Obstacle problems in mathematical physics_ , Elsevier 1987.
  * [50] G. Stampacchia: _Formes bilineaires coercivitives sur les ensembles convexes,_ C.R. Ac. Sci. Paris, 258, (1964), 4413-4416.