Multi-Operator Two-Phase Free-Boundary Problem

• Multi-operator two-phase free-boundary problems are models where distinct PDE operators govern complementary regions coupled by precise interface conditions.
• They incorporate diverse methodologies including elliptic analysis, spectral decomposition, and variational formulations to address interface dynamics.
• Applications span compressible-incompressible flows, geometric measure theory, and tumor growth models, providing insights into stability and symmetry breaking.

Multi-operator two-phase free-boundary problems are problems in which an unknown interface separates two regions governed by different operators, while the interface itself is determined by additional transmission, curvature, kinematic, or measure-theoretic conditions. In current formulations this includes divergence-form elliptic problems with piecewise constant coefficients and mean-curvature overdetermination (Cavallina, 2020), compressible-incompressible flows coupling Navier-Stokes-Korteweg and Navier-Stokes dynamics through phase transition and surface tension (Watanabe, 2018), graph-based reductions in which the free-boundary motion becomes a nonlocal parabolic equation on the interface (Chang-Lara et al., 2018), and blow-up limits for two-phase elliptic measure with frozen coefficients in a multi-operator setting (Goering et al., 4 Sep 2025). The unifying feature is that each phase carries its own operator, while the free boundary supplies the coupling law.

1. Canonical formulations and the meaning of “multi-operator”

A standard two-phase evolution formulation considers a function $U=U(X,t)$, $X\in\mathbb{R}^{d+1}$, with positivity and negativity sets

$$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$

and interface $\Gamma(t)=\partial\Omega^+(t)$. If $L^+=F_1(D^2\cdot,\nabla\cdot)$ and $L^-=F_2(D^2\cdot,\nabla\cdot)$ are translation-invariant, uniformly elliptic operators, the classical evolution is

$$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$

where $G:(0,\infty)^2\to\mathbb{R}$ is Lipschitz and monotone, with $\partial_aG\ge \lambda_0>0$ and $-\partial_bG\ge \lambda_0>0$ (Chang-Lara et al., 2018). In this formulation, the two operators act in complementary phases and the normal velocity is a nonlinear function of the one-sided normal derivatives.

A second prototype arises in compressible-incompressible flow. There the sharp interface

$X\in\mathbb{R}^{d+1}$0

separates a compressible region $X\in\mathbb{R}^{d+1}$1 from an incompressible region $X\in\mathbb{R}^{d+1}$2, and the physical unknowns are $X\in\mathbb{R}^{d+1}$3 together with the moving interface. The compressible phase is governed by a Navier-Stokes-Korteweg operator and the incompressible phase by the Navier-Stokes or Stokes operator, while the interface enforces mass conservation, stress balance, and a generalized Stefan-Gibbs-Thomson law (Watanabe, 2018). In the bounded-domain formulation, the same coupling is described as a free boundary problem of compressible-incompressible two-phase flows with surface tension and phase transition, with maximal $X\in\mathbb{R}^{d+1}$4-regularity and analytic-semigroup structure (Watanabe, 2018).

In the geometric-measure setting, the bulk operators are divergence-form elliptic operators

$X\in\mathbb{R}^{d+1}$5

on complementary NTA domains $X\in\mathbb{R}^{d+1}$6, with symmetric uniformly elliptic coefficients $X\in\mathbb{R}^{d+1}$7. After blow-up at a boundary point $X\in\mathbb{R}^{d+1}$8, one obtains a global free-boundary PDE

$X\in\mathbb{R}^{d+1}$9

together with equality of the conormal derivatives across the free boundary (Goering et al., 4 Sep 2025). This frozen-coefficient problem is a canonical multi-operator two-phase model in the elliptic category.

The phrase “multi-operator” is therefore not restricted to a single PDE class. In one formulation it refers explicitly to coupling of two distinct evolution operators, one for compressible Navier-Stokes-Korteweg in $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$0 and one for incompressible Navier-Stokes in $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$1 (Watanabe, 2018). In another, it refers to complementary elliptic operators with different coefficient fields $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$2 whose interaction is detected through elliptic measure and blow-up analysis (Goering et al., 4 Sep 2025).

2. Elliptic overdetermination with mean curvature

A static two-phase elliptic model is given by the overdetermined problem on a bounded $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$3-domain $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$4, $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$5, containing a $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$6-subdomain $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$7 such that $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$8 is connected. For a contrast parameter $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$9, $\Gamma(t)=\partial\Omega^+(t)$0, one defines

$\Gamma(t)=\partial\Omega^+(t)$1

and seeks $\Gamma(t)=\partial\Omega^+(t)$2 and $\Gamma(t)=\partial\Omega^+(t)$3 such that

$\Gamma(t)=\partial\Omega^+(t)$4

where $\Gamma(t)=\partial\Omega^+(t)$5 is the mean curvature of $\Gamma(t)=\partial\Omega^+(t)$6 and $\Gamma(t)=\partial\Omega^+(t)$7 is the outward normal derivative (Cavallina, 2020). If $\Gamma(t)=\partial\Omega^+(t)$8 are the principal curvatures, then

$\Gamma(t)=\partial\Omega^+(t)$9

and for a ball of radius $L^+=F_1(D^2\cdot,\nabla\cdot)$0, $L^+=F_1(D^2\cdot,\nabla\cdot)$1.

The local analysis is carried out near the trivial configuration $L^+=F_1(D^2\cdot,\nabla\cdot)$2 consisting of two concentric balls of radii $L^+=F_1(D^2\cdot,\nabla\cdot)$3. Perturbations of the inner core and outer boundary are written as

$L^+=F_1(D^2\cdot,\nabla\cdot)$4

with $L^+=F_1(D^2\cdot,\nabla\cdot)$5 and $L^+=F_1(D^2\cdot,\nabla\cdot)$6 in suitable $L^+=F_1(D^2\cdot,\nabla\cdot)$7-spaces of zero-mean functions. The overdetermined condition is encoded by

$L^+=F_1(D^2\cdot,\nabla\cdot)$8

where $L^+=F_1(D^2\cdot,\nabla\cdot)$9 solves the Dirichlet problem in $L^-=F_2(D^2\cdot,\nabla\cdot)$0 with coefficient $L^-=F_2(D^2\cdot,\nabla\cdot)$1 shaped by $L^-=F_2(D^2\cdot,\nabla\cdot)$2, and $L^-=F_2(D^2\cdot,\nabla\cdot)$3 projects onto zero-average functions on $L^-=F_2(D^2\cdot,\nabla\cdot)$4 (Cavallina, 2020). The identity $L^-=F_2(D^2\cdot,\nabla\cdot)$5 is equivalent to the statement that the normal second derivative of $L^-=F_2(D^2\cdot,\nabla\cdot)$6 is constant on $L^-=F_2(D^2\cdot,\nabla\cdot)$7, equivalently $L^-=F_2(D^2\cdot,\nabla\cdot)$8.

The key linearized operator is diagonalized in spherical harmonics. At $L^-=F_2(D^2\cdot,\nabla\cdot)$9,

$$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$0

and the zeros of the explicit eigenvalues $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$1 determine a family of critical contrasts $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$2 and the critical set $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$3 (Cavallina, 2020). When $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$4, the linearized operator is invertible, so the Banach-space Implicit Function Theorem yields a local $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$5 branch

$$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$6

of nontrivial solutions, with $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$7.

The critical case is structurally different. If $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$8, then exactly one spherical-harmonic mode of degree $$L^+U=0 \ \text{in }\Omega^+(t), \qquad L^-U=0 \ \text{in }\Omega^-(t), \qquad V(\Gamma(t))=G(\partial_n^+U,\partial_n^-U) \ \text{on }\Gamma(t),$$9 lies in $G:(0,\infty)^2\to\mathbb{R}$0, and the Crandall-Rabinowitz theorem gives a one-parameter branch of symmetry-breaking solutions

$G:(0,\infty)^2\to\mathbb{R}$1

where $G:(0,\infty)^2\to\mathbb{R}$2 is a nonzero $G:(0,\infty)^2\to\mathbb{R}$3-th spherical harmonic (Cavallina, 2020). Along this branch the inner core remains the unit ball while the outer boundary breaks full rotational symmetry.

The comparison with Serrin-type overdetermination is a central point of interpretation. In the one-phase limit $G:(0,\infty)^2\to\mathbb{R}$4 or $G:(0,\infty)^2\to\mathbb{R}$5, the trivial ball solutions are isolated up to translations. By contrast, the two-phase mean-curvature condition admits nontrivial perturbations and infinitely many nontrivial shapes (Cavallina, 2020). This rules out the common expectation that overdetermined elliptic conditions should enforce spherical symmetry in all two-phase settings.

3. Compressible-incompressible flows with phase transition and surface tension

In the flow-theoretic setting, $G:(0,\infty)^2\to\mathbb{R}$6 and $G:(0,\infty)^2\to\mathbb{R}$7 are separated by a moving interface $G:(0,\infty)^2\to\mathbb{R}$8, with compressible variables $G:(0,\infty)^2\to\mathbb{R}$9 in $\partial_aG\ge \lambda_0>0$0 and incompressible variables $\partial_aG\ge \lambda_0>0$1 in $\partial_aG\ge \lambda_0>0$2. The strong bulk equations are

$\partial_aG\ge \lambda_0>0$3

and

$\partial_aG\ge \lambda_0>0$4

with

$\partial_aG\ge \lambda_0>0$5

in the bounded-domain presentation (Watanabe, 2018). In the general-domain strong-solution framework, the Korteweg stress is written explicitly as

$\partial_aG\ge \lambda_0>0$6

while

$\partial_aG\ge \lambda_0>0$7

in the incompressible phase (Watanabe, 2018).

The interface carries three coupled laws. First, the kinematic or mass-balance condition identifies the normal velocity $\partial_aG\ge \lambda_0>0$8 with the normal traces of the velocities and introduces the phase flux $\partial_aG\ge \lambda_0>0$9: $-\partial_bG\ge \lambda_0>0$0 in the bounded-domain formulation (Watanabe, 2018). Second, the jump of normal momentum is balanced by surface tension,

$-\partial_bG\ge \lambda_0>0$1

with $-\partial_bG\ge \lambda_0>0$2 the $-\partial_bG\ge \lambda_0>0$3-mean curvature (Watanabe, 2018). Third, thermodynamic consistency is imposed through a Gibbs-Thomson or free-energy jump condition. In one formulation,

$-\partial_bG\ge \lambda_0>0$4

(Watanabe, 2018); in another, the phase-transition law is written as

$-\partial_bG\ge \lambda_0>0$5

(Watanabe, 2018).

The geometric reduction to a fixed domain uses a Hanzawa-type transform. In the bounded-domain theory, the reference configuration is a ball $-\partial_bG\ge \lambda_0>0$6, the interface is parametrized by a height function $-\partial_bG\ge \lambda_0>0$7 over $-\partial_bG\ge \lambda_0>0$8, and the map

$-\partial_bG\ge \lambda_0>0$9

flattens $X\in\mathbb{R}^{d+1}$00 to $X\in\mathbb{R}^{d+1}$01 while fixing the outer boundary (Watanabe, 2018). In the general-domain strong-solution theory, one likewise fixes a reference interface $X\in\mathbb{R}^{d+1}$02 and rewrites the quasilinear problem on fixed $X\in\mathbb{R}^{d+1}$03 with interface height $X\in\mathbb{R}^{d+1}$04 (Watanabe, 2018).

The functional-analytic framework is maximal $X\in\mathbb{R}^{d+1}$05 regularity. In the bounded-domain global theory one takes $X\in\mathbb{R}^{d+1}$06, $X\in\mathbb{R}^{d+1}$07, $X\in\mathbb{R}^{d+1}$08, and works with

$X\in\mathbb{R}^{d+1}$09

$X\in\mathbb{R}^{d+1}$10

$X\in\mathbb{R}^{d+1}$11

$X\in\mathbb{R}^{d+1}$12

(Watanabe, 2018). In the strong-solution theorem for general domains, one obtains

$X\in\mathbb{R}^{d+1}$13

for sufficiently small data (Watanabe, 2018).

The main results separate into local-in-time strong solvability in general domains and global solvability near a ball in bounded domains. For every $X\in\mathbb{R}^{d+1}$14, the general-domain problem admits a unique strong solution on $X\in\mathbb{R}^{d+1}$15 in the maximal $X\in\mathbb{R}^{d+1}$16 class provided the initial data are small in their natural norms (Watanabe, 2018). In the bounded-domain setting, if the initial data are sufficiently small and the initial incompressible domain is close to a ball, there exists a unique global solution in the maximal $X\in\mathbb{R}^{d+1}$17-regularity class and the corresponding analytic semigroup is exponentially stable on the infinite time interval (Watanabe, 2018). The analysis exploits analytic semigroups, spectral estimates, uniform $X\in\mathbb{R}^{d+1}$18-bounds for resolvent families, and a nonlinear contraction argument.

4. Interface equations, nonlocal reductions, and variational formulations

One major line of development replaces the moving-boundary problem by an equation posed directly on the interface. If the free boundary is a graph

$X\in\mathbb{R}^{d+1}$19

one solves auxiliary elliptic problems in the upper and lower strips, defines the Dirichlet-to-Neumann-type quantities

$X\in\mathbb{R}^{d+1}$20

and then sets

$X\in\mathbb{R}^{d+1}$21

The free boundary law becomes

$X\in\mathbb{R}^{d+1}$22

in the two-phase case (Chang-Lara et al., 2018). On sets of functions with uniform bounds $X\in\mathbb{R}^{d+1}$23 and $X\in\mathbb{R}^{d+1}$24, the resulting operator admits a min-max integro-differential representation in terms of Lévy operators, and the free-boundary operator has the Global Comparison Property. This yields a viscosity theory for $X\in\mathbb{R}^{d+1}$25, uniqueness by comparison, existence by Perron’s method and barrier functions, and propagation of modulus of continuity.

A different reformulation is variational. In the two-phase quadrature-surface problem one seeks an open-set decomposition

$X\in\mathbb{R}^{d+1}$26

together with a function $X\in\mathbb{R}^{d+1}$27 satisfying

$X\in\mathbb{R}^{d+1}$28

so that for every harmonic $X\in\mathbb{R}^{d+1}$29 continuous up to $X\in\mathbb{R}^{d+1}$30,

$X\in\mathbb{R}^{d+1}$31

(Arakelyan et al., 2016). The corresponding minimization problem is based on the two-phase Bernoulli functional

$X\in\mathbb{R}^{d+1}$32

with $X\in\mathbb{R}^{d+1}$33 and $X\in\mathbb{R}^{d+1}$34 smooth approximations of $X\in\mathbb{R}^{d+1}$35 and $X\in\mathbb{R}^{d+1}$36.

The Euler-Lagrange structure combines bulk equations and free-boundary conditions. The bulk equation is

$X\in\mathbb{R}^{d+1}$37

while domain variation yields the Bernoulli condition $X\in\mathbb{R}^{d+1}$38 on the free boundary. At smooth two-phase points, this is equivalent to

$X\in\mathbb{R}^{d+1}$39

when both phases meet (Arakelyan et al., 2016).

The regularity theory distinguishes one-phase points, two-phase non-branch points, and branch points. One-phase points form a $X\in\mathbb{R}^{d+1}$40-hypersurface. Two-phase non-branch points are real-analytic level sets of the harmonic function $X\in\mathbb{R}^{d+1}$41. Branch points lie on the union of two $X\in\mathbb{R}^{d+1}$42 surfaces meeting tangentially, and the free boundary has locally finite $X\in\mathbb{R}^{d+1}$43-dimensional Hausdorff measure (Arakelyan et al., 2016). In the multi-phase extension, the key geometric statement is that no point can be the common meeting of three or more phases away from the supports of the forcing measures. This is proved through non-degeneracy, local Lipschitz regularity, and a three-phase monotonicity formula.

These two reformulations illuminate different aspects of the same class of problems. The nonlocal parabolic approach emphasizes comparison and viscosity methods on the interface (Chang-Lara et al., 2018), whereas the variational Bernoulli approach emphasizes existence, free-boundary regularity, and junction structure (Arakelyan et al., 2016). A plausible implication is that the operator coupling may be studied either through boundary evolution operators or through energies whose first variation encodes the interface law.

5. Spherical harmonics, mode decomposition, and stability mechanisms

Mode decomposition by spherical harmonics is a recurring tool in multi-operator two-phase problems. In the mean-curvature overdetermined elliptic problem, the linearization in the boundary perturbation variable is diagonal in spherical harmonics, and the eigenvalues $X\in\mathbb{R}^{d+1}$44 determine the critical contrasts $X\in\mathbb{R}^{d+1}$45 at which invertibility fails (Cavallina, 2020). This converts the existence problem into a mode-by-mode spectral analysis and explains why only a single harmonic degree enters the kernel at a simple critical value.

A parallel strategy appears in the free-boundary model of two-phase tumor growth. There the tumor region $X\in\mathbb{R}^{d+1}$46 contains proliferating and quiescent cells, with nutrient concentration $X\in\mathbb{R}^{d+1}$47, proliferating-cell density $X\in\mathbb{R}^{d+1}$48, cell velocity $X\in\mathbb{R}^{d+1}$49, pressure $X\in\mathbb{R}^{d+1}$50, and moving boundary $X\in\mathbb{R}^{d+1}$51. The model combines

$X\in\mathbb{R}^{d+1}$52

with

$X\in\mathbb{R}^{d+1}$53

(Cui, 2013). After linearization around the unique radial stationary solution, perturbations are expanded as

$X\in\mathbb{R}^{d+1}$54

and each mode satisfies a $X\in\mathbb{R}^{d+1}$55 evolution system. Theorem 8.1 then states that there exists $X\in\mathbb{R}^{d+1}$56 and $X\in\mathbb{R}^{d+1}$57 such that for all $X\in\mathbb{R}^{d+1}$58, every non-translational perturbation decays exponentially fast in time; the $X\in\mathbb{R}^{d+1}$59 translation modes are neutral (Cui, 2013).

In the compressible-incompressible flow problem on bounded domains, spectral analysis also isolates the geometrically relevant modes. The generator of the linearized operator has spectrum strictly in $X\in\mathbb{R}^{d+1}$60 on an invariant subspace obtained by excluding zero mass and zero first spherical-harmonic modes of the interface height, and this spectral gap enters the proof of global existence and exponential stability (Watanabe, 2018). Here the decomposition is not merely a technical convenience: it encodes conservation laws, barycenter constraints, and the distinction between decaying and neutral modes.

These results show that multi-operator coupling does not prevent precise spectral resolution. On the contrary, the presence of several operators often makes modal analysis more informative, because each harmonic degree probes a different balance between bulk propagation, interfacial curvature, and transmission. The literature also shows that symmetry breaking and asymptotic stability are not contradictory phenomena: the former occurs at critical loss of invertibility in static elliptic problems (Cavallina, 2020), while the latter occurs in small-data evolutionary regimes after the neutral directions have been factored out (Cui, 2013).

6. Blow-ups, elliptic measure, and geometric structure of the boundary

In the elliptic-measure framework, one starts with an unbounded two-sided NTA domain $X\in\mathbb{R}^{d+1}$61, $X\in\mathbb{R}^{d+1}$62, its complement $X\in\mathbb{R}^{d+1}$63, and common boundary $X\in\mathbb{R}^{d+1}$64. The operators

$X\in\mathbb{R}^{d+1}$65

have symmetric uniformly elliptic coefficients $X\in\mathbb{R}^{d+1}$66 that are $X\in\mathbb{R}^{d+1}$67-quasicontinuous, and give rise to elliptic measures $X\in\mathbb{R}^{d+1}$68 with interior poles (Goering et al., 4 Sep 2025). Writing

$X\in\mathbb{R}^{d+1}$69

the boundary is subdivided into

$X\in\mathbb{R}^{d+1}$70

and one further restricts to a full-measure subset $X\in\mathbb{R}^{d+1}$71 on which $X\in\mathbb{R}^{d+1}$72, $X\in\mathbb{R}^{d+1}$73, and the density of $X\in\mathbb{R}^{d+1}$74 behave well.

The main structural theorem gives a disjoint decomposition

$X\in\mathbb{R}^{d+1}$75

such that for every $X\in\mathbb{R}^{d+1}$76, the tangent measures satisfy

$X\in\mathbb{R}^{d+1}$77

where $X\in\mathbb{R}^{d+1}$78 is the cone of flat measures $X\in\mathbb{R}^{d+1}$79, $X\in\mathbb{R}^{d+1}$80 is the singular part on which $X\in\mathbb{R}^{d+1}$81, and $X\in\mathbb{R}^{d+1}$82 (Goering et al., 4 Sep 2025). A corollary gives a partial answer to Oksendal’s conjecture: the mutual-absolute-continuity set $X\in\mathbb{R}^{d+1}$83 has Hausdorff dimension at most $X\in\mathbb{R}^{d+1}$84.

The route to this result passes through a reduction to a multi-operator two-phase free-boundary problem. At a point $X\in\mathbb{R}^{d+1}$85, one rescales the domains, elliptic measures, and Green functions. Compactness then yields limiting domains $X\in\mathbb{R}^{d+1}$86, limiting measures $X\in\mathbb{R}^{d+1}$87, and limiting harmonic functions $X\in\mathbb{R}^{d+1}$88 with

$X\in\mathbb{R}^{d+1}$89

and $X\in\mathbb{R}^{d+1}$90 (Goering et al., 4 Sep 2025). Writing $X\in\mathbb{R}^{d+1}$91 produces the global free-boundary PDE with frozen coefficients described in Section 1.

The rigidity step classifies globally flat solutions. If a blow-up measure stays uniformly close to the cone of flat measures at all large scales, then the associated global solution is exactly a two-plane solution

$X\in\mathbb{R}^{d+1}$92

so the support of the limiting measure is a hyperplane and the measure is flat (Goering et al., 4 Sep 2025). The proof combines compactness, an $X\in\mathbb{R}^{d+1}$93-monotonicity argument, Caffarelli’s monotonicity-improvement and free-boundary regularity theory, and Preiss’s connectedness lemma for tangent-measure cones.

The geometric consequence is that at almost every point of $X\in\mathbb{R}^{d+1}$94, the boundary admits vanishing Jones $X\in\mathbb{R}^{d+1}$95-number,

$X\in\mathbb{R}^{d+1}$96

hence local planar approximation. Under further hypotheses such as $X\in\mathbb{R}^{d+1}$97 Radon, the good set $X\in\mathbb{R}^{d+1}$98 is $X\in\mathbb{R}^{d+1}$99-rectifiable and $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$00 on $$\Omega^+(t)=\{U(\cdot,t)>0\}, \qquad \Omega^-(t)=\{U(\cdot,t)<0\},$$01 (Goering et al., 4 Sep 2025). This places multi-operator two-phase free-boundary problems in direct contact with tangent-measure theory, rectifiability, and the fine structure of elliptic measure.

Across these formulations, the subject is characterized by a fixed pattern: distinct bulk operators, a free boundary that closes the system, and analytical tools adapted to the interface law. Shape derivatives and implicit-function arguments describe local elliptic branches (Cavallina, 2020); analytic semigroups and maximal regularity control time-dependent two-phase flows (Watanabe, 2018); nonlocal interface equations provide comparison principles (Chang-Lara et al., 2018); variational methods resolve existence and junction geometry (Arakelyan et al., 2016); spherical-harmonic decompositions isolate critical and neutral modes (Cui, 2013); and blow-up theory connects operator coupling to flatness and rectifiability of the boundary (Goering et al., 4 Sep 2025).