Nonlocal Double Phase Operator

• Nonlocal double phase operator is a nonlinear integro-differential operator combining two fractional phases with distinct growth exponents and differentiability orders.
• It employs mixed fractional Sobolev spaces, tail estimates, and dual formulations to bridge weak and viscosity solution frameworks.
• Recent studies reveal self-improving regularity, robust Hölder and gradient estimates, and precise De Giorgi class characterizations under phase-specific conditions.

The nonlocal double phase operator is a nonlinear integro-differential operator that combines two fractional interaction mechanisms with different growth exponents and, in general, different differentiability orders. In its prototype form, it is written as

$$\mathcal L u(x)=\operatorname{P.V.}\int_{\mathbb R^n}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{sp}(x,y)\,dy+\operatorname{P.V.}\int_{\mathbb R^n}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{tq}(x,y)\,dy,$$

or, in the notation of the nonhomogeneous theory,

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$

The first integral is the $p$-phase, the second is the $q$-phase weighted by a modulating coefficient $a(x,y)\ge 0$. This structure is the nonlocal analogue of the classical local double phase operator, and it includes the fractional $p$-Laplacian as the special case $a\equiv 0$ (Ghosh et al., 22 May 2025, Fang et al., 2021).

1. Structural form and phase interaction

The defining feature of the operator is the coexistence of two distinct fractional elliptic phases. One phase has $p$-growth and fractional order $s$; the other has $q$-growth and fractional order $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$0, and is activated or suppressed by the coefficient $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$1. The terminology “double phase” refers precisely to this switching between a $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$2-phase and a $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$3-phase. In the nonlocal setting the switch occurs at the level of pairwise interactions $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$4, rather than pointwise in $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$5 alone (Filippis et al., 2019, Scott et al., 2020).

The kernels are typically assumed symmetric and comparable to fractional kernels. Representative hypotheses are

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$6

or, in the notation of the nonhomogeneous paper,

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$7

Symmetry, continuity away from the diagonal, and, in several papers, translation invariance are also imposed (Ghosh et al., 22 May 2025, Fang et al., 2021).

Two parameter regimes recur in the literature. One is the “lower-order perturbation” regime $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$8, which appears in weak/viscosity equivalence, self-improvement, parabolic Hölder theory, and Harnack theory (Ghosh et al., 22 May 2025, Fang et al., 2021, Giacomoni et al., 2021, Kim, 9 Sep 2025). The other is the harder regime $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$9, where the second phase has the same or higher nonlocal order than the first; in that case Hölder continuity of weak solutions is proved under stronger compatibility conditions such as

$p$0

when $p$1 is Hölder continuous, and gradient Hölder regularity is obtained under

$p$2

when $p$3 is Lipschitz (Byun et al., 2021, Fang et al., 24 Apr 2026).

This split is central to the theory. It shows that “nonlocal double phase operator” is not a single fixed model, but a family of mixed-growth fractional operators whose analytical behavior depends sharply on the relation between the two phases.

2. Variational, weak, and tail-based formulations

The operator is naturally associated with mixed fractional energies of $p$4-type. In the weak theory, one works in fractional Sobolev spaces together with tail spaces that encode the influence of values outside the localization region. A standard requirement is

$p$5

or analogous variants such as $p$6 and $p$7 (Ghosh et al., 22 May 2025, Fang et al., 2021).

For the operator $p$8, the weak formulation is expressed through

$p$9

where $q$0. A weak supersolution satisfies

$q$1

for every nonnegative $q$2, with equality in the solution case (Ghosh et al., 22 May 2025).

The nonhomogeneous theory introduces the nonlocal energy densities

$q$3

so that the right-hand side may depend on $q$4, $q$5, and both phase-dependent nonlocal energies (Ghosh et al., 22 May 2025).

A distinctive analytic feature is the tail. The theory uses quantities such as $q$6 and $q$7, or their $q$8- and $q$9-variants, to measure the contribution of $a(x,y)\ge 0$0 to the operator. This tail dependence is not auxiliary: it enters Caccioppoli estimates, logarithmic lemmas, local boundedness, Hölder estimates, and Harnack inequalities throughout the subject (Ghosh et al., 22 May 2025, Fang et al., 2021, Kim, 9 Sep 2025).

The functional framework also admits rougher generalizations. One paper replaces the prototype kernels by measurable kernels $a(x,y)\ge 0$1 satisfying two-sided ellipticity bounds and nonlinearities $a(x,y)\ge 0$2 satisfying

$a(x,y)\ge 0$3

and proves the same self-improving conclusions for this broader class (Scott et al., 2020).

3. Viscosity solutions and equivalence with weak solutions

The viscosity framework adapts the standard touching-test paradigm to a nonlocal mixed-growth operator. In the homogeneous setting, a viscosity solution is tested with a $a(x,y)\ge 0$4 function patched to the original solution outside a contact ball; in the nonhomogeneous setting the test inequality must also evaluate the phase-dependent quantities $a(x,y)\ge 0$5 and $a(x,y)\ge 0$6 on the test function (Fang et al., 2021, Ghosh et al., 22 May 2025).

For the nonhomogeneous operator, a viscosity supersolution is lower semicontinuous, satisfies the appropriate tail integrability, and obeys

$a(x,y)\ge 0$7

whenever $a(x,y)\ge 0$8 touches $a(x,y)\ge 0$9 from below at $p$0, under one of the regularity regimes required in the paper: either $p$1, or $p$2, or $p$3 with $p$4 an isolated critical point and $p$5 for $p$6 (Ghosh et al., 22 May 2025).

A foundational pointwise compatibility statement was proved for bounded viscosity solutions under the natural energy finiteness assumption

$p$7

Under the structural assumptions of that paper, if $p$8 in viscosity sense and a $p$9 function touches from above, then $a\equiv 0$0 exists pointwise and satisfies the same inequality at the contact point (Filippis et al., 2019).

The first systematic weak/viscosity bridge for the nonlocal double phase equation was established in the homogeneous setting. Bounded weak supersolutions were shown to be viscosity supersolutions, using local Hölder continuity, a comparison principle, continuity of the operator on smooth functions, and the fact that affine functions are annihilated under the symmetry and translation-invariance assumptions (Fang et al., 2021).

The nonhomogeneous theory then established a two-way equivalence. If $a\equiv 0$1 is a continuous weak supersolution, a suitable comparison principle holds, and $a\equiv 0$2 is continuous in $a\equiv 0$3 and $a\equiv 0$4 and Lipschitz in the nonlocal derivative variables, then $a\equiv 0$5 is a viscosity supersolution. Conversely, if $a\equiv 0$6 is a bounded viscosity supersolution and $a\equiv 0$7 is uniformly continuous in $a\equiv 0$8, Lipschitz in the nonlocal derivative variables, nonincreasing in $a\equiv 0$9, and satisfies the stated growth condition, then $p$0 is a weak supersolution. The second implication uses infimal convolution $p$1, an approximating source $p$2, convergence of the nonlocal energy terms, and passage to the limit in the weak formulation (Ghosh et al., 22 May 2025).

A persistent misconception is that weak and viscosity notions are automatically interchangeable for this class of operators. The published theory shows instead that equivalence is a theorem with nontrivial hypotheses, not a formal identity. In particular, one paper explicitly noted that the reverse implication viscosity $p$3 weak remained open in its setting, and a later paper resolved it for a nonhomogeneous class under additional structural conditions on $p$4 (Fang et al., 2021, Ghosh et al., 22 May 2025).

4. Regularity theory

The regularity theory begins with Hölder continuity. The first Hölder result for bounded viscosity solutions proved that if

$p$5

then $p$6 for some $p$7, assuming only boundedness of the modulating coefficient

$p$8

and suitable exponent restrictions. The proof follows a nonlocal oscillation-decay scheme in the style of Krylov–Safonov rather than a De Giorgi–Nash–Moser energy argument (Filippis et al., 2019).

The weak-solution theory in the $p$9 regime was then developed by De Giorgi–Nash–Moser methods. Under the structural assumptions of the homogeneous equation, bounded weak solutions are locally Hölder continuous, and the oscillation satisfies

$s$0

for $s$1 and $s$2 (Fang et al., 2021).

A different phenomenon appears in the higher-order regime $s$3. For weak solutions in that setting, local boundedness holds under the growth restriction

$s$4

or without a restriction on $s$5 when $s$6. Local Hölder continuity is obtained if, in addition, the coefficient is Hölder continuous and

$s$7

That paper presents these assumptions as the nonlocal analogues of the sharp local double phase conditions (Byun et al., 2021).

The parabolic theory extends the picture. For

$s$8

higher space-time Hölder continuity of weak solutions is proved under $s$9, $q$0, local translation invariance of $q$1, boundedness of $q$2, and the corresponding tail assumptions. A stationary counterpart gives local Hölder continuity under diagonal continuity of $q$3, and a global boundary result yields

$q$4

when $q$5 is $q$6, the right-hand side is bounded, and $q$7 with $q$8 (Giacomoni et al., 2021).

The current highest regularity result concerns gradients. Interior $q$9 regularity of viscosity solutions has been proved for the degenerate case, where $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$00 may vanish, assuming $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$01 is symmetric, translation invariant, bounded, and Lipschitz, and that

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$02

That work explicitly states that it resolves the higher regularity issue raised by De Filippis–Palatucci (Fang et al., 24 Apr 2026).

A common oversimplification is to say that boundedness of $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$03 is either always sufficient or always insufficient for regularity. The published results show a more nuanced picture. Bounded measurable $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$04 suffices for Hölder continuity of bounded viscosity solutions in one regime (Filippis et al., 2019), whereas Hölder continuity of $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$05 is required for Hölder regularity of weak solutions in the $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$06 regime (Byun et al., 2021), and Lipschitz continuity of $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$07 is used for $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$08 regularity (Fang et al., 24 Apr 2026). This suggests that the threshold on the coefficient depends both on the phase interaction and on the level of regularity being sought.

5. Self-improvement, De Giorgi classes, and Harnack theory

Beyond baseline Hölder continuity, the operator exhibits self-improving regularity. For bounded weak solutions under

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$09

and with $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$10, the energy density improves in integrability and the solution gains differentiability and integrability: $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$11 for some $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$12. The proof uses a dual-pair reformulation $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$13, a reverse Hölder inequality on diagonal sets, and a fractional Gehring lemma adapted from Kuusi–Mingione–Sire (Scott et al., 2020).

A related self-improving theory has been proved for nonlocal double phase equations with VMO coefficients and solution-dependent principal coefficient $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$14. Under $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$15, weak solutions satisfy a local higher differentiability and integrability property,

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$16

and under stronger assumptions they satisfy

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$17

where

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$18

The mechanism there is perturbative freezing of VMO coefficients together with Campanato iteration (Byun et al., 2023).

The nonhomogeneous theory also develops a double phase De Giorgi class $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$19 for the simpler equation

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$20

These classes encode Caccioppoli-type inequalities simultaneously for the $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$21-fractional seminorm, the weighted $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$22-fractional seminorm, and the nonlocal tails. Weak solutions in the relevant Sobolev and tail spaces are then locally bounded in $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$23 (Ghosh et al., 22 May 2025).

Potential-theoretic estimates have now reached Harnack theory. For weak solutions with positive bounded symmetric modulating coefficient and no Hölder condition on $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$24, local boundedness, weak Harnack inequalities, and full Harnack inequalities have been proved under $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$25. For nonnegative weak solutions,

$$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$26

and in the globally nonnegative case the tail term disappears (Kim, 9 Sep 2025).

Taken together, these results show that the nonlocal double phase operator supports a regularity theory extending from boundedness and Hölder continuity to higher Sobolev gain, De Giorgi classes, and Harnack estimates, with the exact conclusions controlled by the relation between the two phases and the regularity of the coefficient.

6. Extensions, variants, and boundaries of the notion

Several later developments extend the framework without changing its core double phase principle. One direction mixes local and nonlocal terms. A parabolic model couples a local $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$27-growth diffusion $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$28 with a nonlocal $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$29-Laplace type operator, and proves local boundedness of weak solutions by De Giorgi–Nash–Moser iteration; under $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$30 and $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$31, it also proves semicontinuity and pointwise behavior of supersolutions (Shang et al., 2023). Another work studies mixed local/nonlocal double phase equations in a variational setting, including a model with a local $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$32-Laplacian coupled to a fractional double phase operator and proves existence, constant-sign solutions, and a least energy sign-changing solution (Arora et al., 1 Mar 2026).

A second direction introduces more elaborate Musielak-type structures. One paper studies a variable-order Neumann and Robin problem built from the sum of two nonlocal fractional Musielak operators with orders $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$33 and $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$34, develops the associated fractional Musielak–Sobolev spaces, and proves existence of weak solutions by Ekeland’s variational principle and direct minimization (Srati, 2024).

At the same time, the literature contains nearby notions that should not be conflated with the nonlocal double phase operator proper. A local anisotropic framework with exponents depending on $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$35 or on $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$36 develops a new local double phase operator in Musielak–Orlicz spaces, but it explicitly does not involve a nonlocal fractional operator (Bahrouni et al., 2024). Likewise, non-isothermal phase-field models with a double-obstacle potential use a bounded nonlocal convolution operator $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$37, yet the “two-phase” aspect there comes from the obstacle potential rather than from a two-growth integro-differential operator (Burkovska, 2023).

The modern theory therefore treats the nonlocal double phase operator as a distinct object: a mixed fractional $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$38- and $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$39-growth operator, usually with a modulating coefficient $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$40, whose analysis combines nonlocal tails, mixed-growth energies, De Giorgi or oscillation methods, comparison principles, and, increasingly, refined perturbative and linearization techniques. The cumulative results show that weak and viscosity theories can be made compatible, that regularity depends sharply on the phase balance, and that the operator supports a broad analytic program ranging from self-improving estimates to Harnack inequalities and interior $$L_a u(x)=2\,\operatorname{P.V.}\int_{\mathbb R^N}|u(x)-u(y)|^{p-2}(u(x)-u(y))K_{s,p}(x,y)\,dy+2\,\operatorname{P.V.}\int_{\mathbb R^N}a(x,y)|u(x)-u(y)|^{q-2}(u(x)-u(y))K_{t,q}(x,y)\,dy.$$41 theory (Ghosh et al., 22 May 2025, Fang et al., 24 Apr 2026).