Efficient Differentiable Contact Model with Long-range Influence
Authors: XIAOHAN YE
Kui Wu
Taku Komura
zherong pan
Venue: ICLR 2026 Poster (2026)
Keywords: Differentiable Simulation
Policy Search
Motion Planning
Primary area: applications to robotics, autonomy, planning

Abstract: With the maturation of differentiable physics, its role in various downstream applications—such as model-predictive control, robotic design optimization, and neural PDE solvers—has become increasingly important. However, the derivative information provided by differentiable simulators can exhibit abrupt changes or vanish altogether, impeding the convergence of gradient-based optimizers. In this work, we demonstrate that such erratic gradient behavior is closely tied to the design of contact models. We further introduce a set of properties that a contact model must satisfy to ensure well-behaved gradient information. Lastly, we present a practical contact model for differentiable rigid-body simulators that satisfies all of these properties while maintaining computational efficiency. Our experiments show that, even from simple initializations, our contact model can discover complex, contact-rich control signals, enabling the successful execution of a range of downstream locomotion and manipulation tasks.

Published as a conference paper at ICLR 2026
EFFICIENT DIFFERENTIABLE CONTACT MODEL WITH LONG-RANGE INFLUENCE
Xiaohan Ye', Kui Wu², Taku Komura' Zherong Pan
The University of Hong Kong LLIHHTSPEED Meta
ABSTRACT
With the maturation differentiable physics, its role in various downstream ap- plications-sur preel-predictive control, robotic design optimization, and neural PDE solvers-has become increasingly important. However, the deriva- tive information provided by differentiable simulators exhibit abrupt changes or vanish altogether, impeding the convergence of mradeent-based optimizers. this work, we demonstrate that such erratic gradient behavior is closely tied the design of contact models. We further introduce a set of properties that a contact model must satisfy to ensure bell-behaved gradient information. Lastly, we present practical contact model for differentiable rigid-body simulators that satisfies all of these properties while maintaining computational efficiency. Our experiments show that, even from simple initializations, our contact model can discover complex, contact-rich control signals, enabling the successful execution of range of downstream locomotion and manipulation tasks.
INTRODUCTION
Recent advancements in differentiable physical models (Werling et al. 2021 Huang et al. 2024) have unlocked range of downstream applications, including model- reinforcement learn- ing Xu et al. (2022). shooting- based controller optimization (Amos et al. 2018). and robot co- design (Xu al.) 2021). State-of-the-an models now extend to various material types, encompassing both rigid and deformable bodies, while offering first-order gradient information. A notable advan- tage of differentiable physics is ability to automatically discover contact-rich motions from trivial initializations (Mordatch et al. 2012 Pang et al. 2023). Achieving this, however, requires an ideal contact model that strikes balance between accurately approximating physical contact mechanisms and providing meaningful gradient information. Over the years, significant efforts have been made to optimize this balance (Werling et al. 2021 Huang et al. 2024 Le Cleac' h et al. 2023). Despite their significant progress, recent systematic analyses (Antonova al. 2023 Suh al. 2022a b)
have highlighted several pitfalls in the gradient information provided by differentiable physics sys-
tems. While analytic gradients are beneficial across much of the objective landscape, they can hibit rugged behavior when optimizing over non-smooth interactions and may vanish in nearly flat regions. Consequently, optimizers are often prone becoming trapped in undesirable local optima.
To address these challenges, existing techniques (Antonova et al. 2023 Li et al. 2022a) employ global optimization algorithms, such as Bayesian optimization cape local optima. While we agree that global search methods are crucial for complementing local gradient-based optimization, we argue that many of the issues with poor gradient information can be mitigated by improving the contact models within existing differentiable physics frameworks. Specifically, both rugged and vanishing gradients stem from the contact model itself. When two rigid objects come into contact, the abrupt introduction of contact forces results in rugged gradients. Conversely, the absence of contact, the lack of direct interaction leads to vanishing gradients.
optimal transportation,
In this paper, we make both theoretical and practical contributions, both aimed at enhancing the gradient landscape of differentiable rigid body simulators. Theoretically, we introduce in Section B set of properties that define well-behaved contact model. These properties ensure that the tact model supports differentiation and can provably prevent inter- pentrreenttraionn generate physically plausible contact forces, and provide non-vanishing gradients even when objects are arbitrarily far

---

R3V The vertices Published as a conference paper at ICLR 2026
apart. The last property leads to long-range influence, allowing us to discover novel contact-rich motions, even from a trivial initialization where objects are distant from each other. Practically, we present computationally tractable contact model in Section | which satisfies all these properties as proved in our Appendix A..1 and is applicable to arbitrary articulated bodies represented using significantly improve the computational efficiency for
triangle meshes. Further in Section evaluating the contact model using a Bounding Sphere Hierarchy (BSH). We have incorporated our contact model into full-featured rigid body simulator and experimented on a row of robotic manip- ulation and locomotion tasks. Our results show that our method can discover complex, contact-rich rich control signals from trivial initialization, while previous models can get optimizers stuck trivial local minima or suffer from slow convergence.
RELATED WORK
The idea of differentiable physics originates from the pioneering work Todorov (2011). which is then built into the MuJoCo simulator (Tassa et al. 2012),, where gradients are computed via costly finite difference schemes. Differentiable simulators are then proposed to use more efficient analytic gradient information. Early works apply this idea to rigid body simulators (de Avila Belbute-Peres et al. 2018) and reduced-order deformable body simulators (Pan & Manocha 2018). The idea is then adopted in other simulator models (Newbury et al., 2024) for elastic and plastic deformations, articulated bodies, and fluid bodies, to name just few. Since their invention, differentiable simu- lators have found many applications in computer graphics, robotics, and machine learning. Early works along this line use differentiable simulators perform model-predictive control (Tassa al.) 2012) and guide deep policy search (Levine & Koltun (013). Diffcrentiable simulators can be combined with appearance models to perform state estimation (Ma al., 2022) and system iden- tification (Le Lidec et al. 2021). In computer graphics. animators use differentiable simulators to inversely design initial and boundary conditions (Li et al.] [023] Stuyck & Chen 2023) Du et al.] 2021). They can also provide gradient information for physics-informed machine learning such as neural PDE solvers (Heiden et al. 202T) and neural motion planning (Toussaint et al. 2019).
Certain substeps in a simulation procedure are inherently non-differentiable, of which the most important substep is contact handling. It is known that gradient information is lost for collision detection with thin-shell-like objects (Harmon al. 2009 Li et al. 2022a bb) and the sudden change of contact forces in collision responses incur non-smoothness, which requires manually choosing a specific direction in the Clark subdifferential (Werling al. (021). These factors
compromise
the gradient information, hindering downstream optimizers performance. In [Suh et al. (2022a), authors propose mixed gradient-free and gradient- based optimizer boost the performance of policy search. On parallel front, [i et al. (2020) showed that the contact model can be reformulated to prevent gradient vanishing for lhin-shell-like objects. Differentiable simulator (Huang et al., 2024) built on top of this technique exhibits better robustness. However, we show that even with this technique, optimizers can still suffer from vanishing gradients.
3 DIFFERENTIABLE PHYSICS WITH WELL-BEHAVED CONTACT MODEL
In this section, we first formulate the problem of a differentiable physical model, and then formalize the properties of a well-behaved contact model that provides useful gradient information.
3.1 DIFFERENIABLE PHYSICS MODEL
Throughout the paper, we consider articulated bodies geometrically represented using triangle meshes. This is the representation adopted by a majority of differentiable contact models (Wer- ling et al. 2021 Huang et al., 2024| Xu et al. 2021). Formally, we assume the configuration of an articulated body is represented using set of V vertices located at Ii.--.V R3, and we x without subscript to denote the concatenated vertex vector connected to form a set of T triangles f1,-,,7, with each cornering three vertices and defined as t; = ((),, ((3) = {1...., V). The configuration space x é C c RaV can be divided into a penetrating set Cutss (ar 3 # ly CH( ) CH(A)A()))) $ Ø), where CH is the closed
2

---

Published as a conference paper at ICLR 2026
convex hull of a set of vertices, and a penetration-free set Ciee C/Cobs. The goal of contact model is to impose contact forces on the rigid object to ensure that x E Cree.
A discrete-time physical model can be cast a time transition function x²+1 f(x', _²-1, ,8),, where we use superscript to denote the time index, is the kinematic state of a body at the tth time instance, and öt is the timestep size. The concatenation of two kinematic states ==²²) composes a dynamic state of the body, with velocity approximated as (x'-x" / It is well-known (Marsden
& West 2001 Gast al. 2015) that the transition function can be cast as a numerical optimization:
argmin x!-1 C(x x ,ô))
which leads to stable performance of modern differentiable position based dynamics, such as Huang al. (2024). In these position-bast models, the Lagrangian function L contains various terms that model different behaviors. Specifically, we define:
C(x x ō() I(x f+1 x², x 5t) + +*P((++*)) D(x*-) ",",,,,,',00 ,t",
where the term models inertial acceleration, P models the contact potential weighted by a param- eter J4, and DD models frictional damping. We refer readers to Huang al. (2024) for more details of other terms and we focus on the term P in this paper. In optimization-
simulators, the contact mechanics are mainly determined by the potential P, whose proper definition has been discussed, e.g., in Fisher & Lin 20001); Harmon et al. (2009): Li et al. (020)) ::000:0000 Ye et al. (2025).
3.2 PROPERTIES OF CONTACT POTENTIAL
We propose four aspects of indispensable properties of a well-behaved contact potential. Barrier potential: It is well-known in the theory of continuous collision handling (Brochu et al. 2012) that penetrating state x' e Cobs corresponds to non-smooth landscape of the contact potential. Therefore, well-behaved contact potential should provably prevent any penetrations by ensuring that x' Cieee for any t. This property is first proposed and achieved in Harmon et al.] (2009). where authors designed P be positive layered potential function that is infinite when x Cohe- As result, numerical optimizer with globalization techniques, such as line search and trust-region, can ensure the monotonic decrease of the Lagrangian C, leading to the finiteness of and thus x'*1 Crice Li itt al. (2020) further elucidates that a potential P should act as a bar-barrier function in interior point optimization, which is formalized as our first property below:
Property 3.1 (Barrier-Form). ((x) ) 0 is continuous for any E Cnd P(x)) oo iff x E Cobs-
Note that Barrier-Form does not describe the exact contact mechanics, because an exact contact model can only impose contact forces on bodies when they are exactly touching, i.e., x E ÖCotss but our contact potential induces the generalized contact force force-UPP /0t even when r Cubs- To mitigate this issue, [i et al. (2020) proposes to iteratively approximate the true contact mechanics by tuning the coefficient/ μs Indeed, we can easily verify that lim μuP converges to the indicator function that equals to oo if x E eCoby and otherwise. Smoothness: In order to utilize the primal log-barrics barrier method for computing z²11 by solving optimization would require P to be at least differentiable.
To this end. [ii et al. (2020) proposed a differentiable surrogate of the triangle-triangle distance function. Unfortunately. although differentiability is enough for solving Equation it is not enough for providing reliable gradient information. Indeed, the gradient of a numerical optimization takes the following form by the implicit function theorem:
On I+1
JC
8C Ox a a 1'
0( x
0.1+
whose proper evaluation requires P to be twice-differentiaale, which not satisfied in Li et al. 2020); Huang al. (2024) as shown in our Appendix A.44 leading to our second property:
Property 3.2 (Smoothness). P is twice differentiable at é Cjree-
Unfortunately, being numerically well-defined does not guarantee that the gradient information can effectively guide the optimizer to find meaningful solutions for downstream applications. To this end, we introduce two other properties that ensure the contact model is non-prehensile and non- vanishing. Non-prehensile: We know that passive contact can only impose unilateral pushing
3

---

Published as a conference paper at ICLR 2026
forces between a pair of contacting objects, instead of pulling objects together. Formally, we in- troduce a sufficient condition to ensure non non-prehenn -oneeeeeeseeee forces. Let us define two index subsets of well-separated vertices I n J 2 and I u J _ {1, V), such that the convex hull of these sets of vertices are non-overlapping. i.e. CH( n CH( = 2. A contact potential between these two subsets can be defined as a pair-wise contact term pluJ (IteI,Fjer) or pTuJ for short. Note that we can also establish [Bariie-Form) and Smoothness for a pair-wise contact term pruj, To this end, we define CTuJ = {z & C|3t, # t2 44 Uutt, S I U J, CC( n CH(Tyet,) * Ø and Barrier-Form and Smoothness for pairwise contact
"bbs
Ch = C//CTs terms:
ohs and we have the following
Definition 3.1( ( Saarier-Formpee
Smoothness for pairwise contact p.uJ pertains Barrier- Form if pDuJ 0 for any E C and paud iff E Cols plud pertains Smoothness if is twice differentiable at Cree
pDud induces the following contact force on any JieI T.jeJ:
FTU
-Tuu dI
(ZIEI,Ize)) feg
-pTuj (IIII, Tza7).
dzzer
To allow only non-prehensile forces, we require that each fuu is pointing from CH(I)e.) to CH( and vice versa. Formally, we define the set of non-zero vectors pointing from CH((T to CH(Ite)) as F--I = {a(a- -b]la Osae CH(TæI) AbE CH(z)er)); and require that:
VieI:f. f.TJ Fy-1 and Vje J SJUJ EFI-J-
(2)
Non-vanishing: Our final property is of paramount importance and ensures that differentiable
simulator provides non-zero gradient information at arbitrary contiguration. This is ensured by definition of the non-prehensile force set Fz-I- Indeed, since CH(rje)) ) and CH( disjoint, closed convex sets, for any (x)ab)= e FI-I, we have # b and a 0. leading to L. for all i ε Z. We further ensure that the contact forces between every pair of geometric primitives (triangles) are taken into consideration. Put together, we can ensure both properties by requiring that P is summation of pairwise contact terms pzud between well-separated vertex clusters:
Property 3.3 -Non-prchensie & non-vanishing). At every x E Crree. we can define a finite family of set pairs A(z) ((I, (I,T)IIn have P E(I.JreA(E) plud such that every term Piud satisfy Equation 2 and for every pair of triangles (t,,tz) different rigid bodies, have t, Ut S I for at least one (I,J) A(x).
2 AIIIUT = VV (1) CH(I n CH(()ee) 2 We
This is an important property that allows the gradient information to be provided for arbitrarily dis- tant objects. In many applications. such gradient information can help a local optimizer discover contact-rich motions from trivial initial guesses. Regretfully, we are not aware of any contact model that pertains Form) [Smoothness] and NNon-prehensile & Non-vanishingl at the same time. We summarize the failure cases for various properties in Figure 2 and compare the property com- pleteness in Table In Figure [ we illustrate the main idea behind our contact model in Section that satisfies Non-prehensile & Non-vanishing
Figure 1: A 2D illustration of contact force between pair of hexagons satisfying Non- prehensile & Non- vanissinii Each hexagons have line segments 2D (resp. triangles in 3D). Left: We compute the exact contact force (arrow) between each pair of nearby line segments. Mid- dle: For faraway pairs of line segments, computing exact contact forces would involve too many segment pairs, e.g. 4 pairs of forces between 2 edges on each hexagon. Right: Instead, we group faraway segments and approximate the contact forces between centers of bounding circles in 2D (resp. bounding spheres 3D).
4

---

Published as a conference paper at ICLR 2026
B Form
Smoothness
Non-prehensile
Non-vanishing
Fullure Case
Figure 2: We illustrate failure cases for various properties, where we assume two brown boxes are separated by distance z and plot the contact potential (x))) that pertains the property in green and fails the property in red. Our BBarrie-Form requires P to tend to infinity as x + 0* [Smoothness] requires P to have well-defined second-order derivatives. [oon-prehensile to always push the two boxes apart. Non-vanishing requires P and thus the contact force to be non non-vanishing for arbitrarily large x.
the contact force
Formulation
Barrier-For Form
Smoothness
Non-prehensile
Non-vanishing
Turpin al. 2022
Sehwarke al.
2025
Werling al.
2021
Xu al.
2022
X
X
Fisher & Lin 2001
Guendelmun et ml. 2003
X
X
X
Harmon ct al
2009
Li et al. 2020
X
/
X
Ye ct al.
2025
/
/
X
Ours
Table I: Comparison of property completeness. Contact models based on complementary condi- tions Werling et al. (202T); Xu et al. (022) or soft penalty functions Fisher & Lin (20T); Guendel- man et al. (2003) cannot guarantee intersection free or sufficient smoothness. Contact models based on the log-barrier functions Harmon et al. (2009); Li et al. (2020) only have first-order derivatives, which does not support differentiation using the inverse function theorem. Finally, prior contact models have vanishing gradient when the distance is larger than small margin. Although, Turpin eet al. (2022); (Schwarke et al. (2025) provides non-vanishing gradients, they still cannot guarantee intersection-free and smoothness, limited by penalty models.
4 WELL-BEHAVED CONTACT POTENTIAL
In this section, we propose practical and well-behaved contact potential. We start by showing that slightly modifying an existing contact potential (Liang et al.] 2024] Ye et al. (025)) makes well-behaved, but such a potential is slow to compute. We then improve its computational efficacy in Section by borrowing ideas from the well-known hierarchical algorithm (Barnes & Hut 1986)
for N-body simulation.
Barrier-Form requires that our contact potential acts as a primal barrier function. However, ex- ssting barrerr potential functions (Harmon et al. 2009 Li et al. 2020) is derived from a modified triangle-triangle distance function, denoted as di (((())),, (Akket,, ). which is then assembled to form the following contact potential: P En,,,, P(d(zn(s))), : ((())et, )). with P being some locally supported barrier function. Unfortunately. this potential is at most first-order differentiable and vio- lates [Smoothness and the local support of P violates Non-vanishing Instead, we propose to adopt the more general contact potential (Liang et al. 2024 Ye et al. 2025) between pair of convex hulls. These methods use the separating hyperplane theorem to prevent two convex hulls from colliding by inserting separating hyperplane between them. This hyperplane is then modeled as a physical object with zero-mass, which serves as auxiliary variables to formulate the contact potential. It has been shown that such general contact potential is globally twice differentiable, which serves good starting point for our derivation. Since triangle is convex by nature, the more general contact potential can be adopted to serve as pl.ut, Specifically, given a pair of triangles t; and ty. since
the two triangles are both convex sets, we can define a separating plane Pij = (n. dz) ) e IR4 between them if the two sets are disjoint, with niy and dej being the normal and offset, such that: (Fa(k))e,,)ni)) ddi > O and (an(())e,, niy) + dij < O. As a result, we introduce the following potential
5

---

Published as a conference paper at ICLR 2026
via nested optimization:
plut,
=min Lig (Piy Fi(k))t, Pjj
Lip(pis Pi(k)et (((H))))))))))
Σ 12 (1- In R-1 ((Ta))) Big) + diy)" R-1 ((z)(1)) -nis) -ds)"
1
3 Σ
where we purposefully introduce a constant coefficient 12 in the first term so that our follow-up derivations take a simpler form, and other positive coefficients can be used well. As a main
point of departure from the original formulation in Liang et al. (2024); Ye et al. (2025). we do not
maxxx, 0) that has global support on R+ to prevent vanishing gradient. Note that naj. daj.dey are computed by minimizing Lis and thus nij not normalized. is easy verify that the objective function defined in pAyut is a strictly convex function with unique minimizer, so that phuee is well-defined function. With the pair-wise potential defined, we can assemble them and define:
locally supported log-barrier function, but define the potential function 1/(0)
P= Pliut (I GA4
(3)
where the summation is taken triangle pairs different rigid bodies. We now show that the so-defined contact potential pertains all our desired properties.
Lemma 4.1. Each pair-wise potential plut, in Equation 3 pertains Barrier-Forn Smoothness
and satisfies Equation 2 that the potential P pertains Barrier- Form Smoothness and Non- prehensile & Non-vanishing
At this point, we have shown that Equation is a well-behaved contact potential function. Re- markably, this function is computationally practical. Indeed, each term pl.ut involves a small 4D- optimization problem with a strictly convex objective function, which can be solved efficiently using Newton's method to evaluate Pij. The first and second derivatives of plutt can then be computed using the inverse function theorem. However, brute force computation of the potential function P is not efficient, since it involves terms that account for the contact potential between each pair of disjoint triangles which increases in the square order of triangles.
5 EFFICIENT CONTACT POTENTIAL EVALUATION
The computational challenge of evaluating Equation 3 lies in accounting for the contact potentials between all pairs of disjoint triangles. This scenario closely parallels the N-body simulation prob- lem, where the forces between all pairs of particles must be computed. Instead of performing a brute-force summation with a computational cost of O(N2), efficient algorithms such as the tree code (Barnes & Hut 1986) and the fast multipole expansion (Greengard & Rokhlin 1987) achieve costs of ONN log N 1) and ONN), respectively. These methods rely on the multipole expansion to separate the influences of source particles from those of target particles. However, several factors make these algorithms unsuitable for our case. First, the multipole expansions for our contact poten- tial pl,ut remain undefined. Second, even if such expansions could be derived, the abrupt transition between the exact potential and its multipole approximation could introduce discontinuities, thereby violating [Smoothness] Inspired by the fast multipole method (Greengard & Rokhlin 1987). propose instead modified potential that is also well-behaved and can be evaluated hierarchically. Our main idea is to smoothly transit from the exact potential function plutt simplified functions that can be hierarchically evaluated.
5.1 SMOOTH TRANSITION BETWEEN POTENTIALS
Let us consider the pairwise potential between index set and J. For an index set I, we define its bounding sphere be centered at JI Σiez wx/[Z] with radius Rr 2 maxar |Iz zz). Suppose there two versions of the potential denoted as PJUJ and pIuJ we can smoothly blend the two functions when the distance between 2y and Iı grows from di to d2 with di d2 as illustrated in Figure yielding the following blending potential:
DIuJ d,2
(1 th (x) (Z) IuU Od -o (x)
(4)

---

Published as a conference paper at ICLR 2026
where we define the interpolation function as:
ody-dd $(([=z= xz -d1)/(d2 d:)) and ((d) max(min(6dd - 15d" 10d 1).
Similar to the tree code algorithm (Barnes & Hut 1986), our goal of blending is to gradually replace exact potential functions with faster-toccompute approximations. Such blending should not happen when two sets of vertices are too close to each other. practice, we only allow blending when RI Rg d1 where Rz, Rg are the radii of the bounding spheres of I.J. We first show the well-behaved nature of potential functions is invariant to blending:
Lemma 5.1. Taking the following assumptions: i) RI + Rg dd < d2; ii) Piy pertains Barrier-
Form
Smoothness and satisfies Equation iii) 0 s ddg S PPh when 1 - xJ di; iv)
d2
has
Smoothness and satisfies Equation then P has the same properties as Phu..
Lemma 5.1 can be immediately used blend our potential function phur, with much simpler, closed-form function. Consider moving all three vertices of I to the center point 2y, (Ia(1) (2) (i(s)) 3/3 and similarly moving ły to Zt,, then the potential plut, takes the following (centered) form after some basic algebraic manipulation:
Plu min
12 (1
12 1+
(5)
(12 may) +g))"
Rsd ((re, -my)) d,,''
1/2
We now ready apply Lemma 5.1 blend pl.tt and puu in well-behaved manner:
Corollary 5.2. If we define Plitt, = plutt and plul, plull, then Pluu, di pertains Barrier-
Form
Smoothness and satisfies Equation
Intuitively, Corollary 5.2 allows us to the exact potential plut, when two triangles very close each other, while switching to simpler, closed-form potential Pi: when the centers of two triangles are well-separated by some distance d2.
5.2 HIERARCHICAL POTENTIAL BLENDING
In the previous section, we employed blend- ing techniques to smoothly replace the costly contact potential plut, with a more computa- tionally efficient closed-form potential p However, since this blending is applied only to individual pairs of triangles, the approach still requires summing over ( terms. In this section, we fully unlock the potential of blend- ing by hierarchically merging triangles to con- struct BSH (Agarwal et al. 2004, Bradshaw & O'Sulliian 2004) for each rigid body, and then smoothly replace the contact potential of each sphere with a single term. Specifically, we adopt layered hierarchy, where each sphere tightly encapsulates all the bounding spheres of its two children. Another widely used option is to use a wrapped hierarchy, where each sphere tightly encapsulates the actual geometry. Although wrapped hierarchy achieves tighter bound, the layered hierarchy is required in our method to ensure well-behaved properties. Our BSH is defined below:
Illustration of our BSH-based contact potential. When two triangles are nearby, we use the exact potential based on separating plane Pij (left). When the center of bounding sphere is sep- arated by at least (Re, Re, (1, c), we use the centered potential in Equation [ (middle). These two cases are combined by smooth blending. The centered potential can be calculated hierarchically for clusters of triangles (right).
Figure 3:
Definition 5.3. A BSH is binary tree, where each node contains an index subset of vertices I g [.,"..,N that is the union of the two subsets of its left and right children, denoted as I I UIr. The radius Ri is the smallest radius encapsulating the bounding spheres of two children. Each leaf node stores a single triangle t.. Further, each node's 's sphere is centered at Iyg with radius Rz.
Throughout the paper, we use the associated index subset to refer to a BSH node. There are many ways to practically construct our BSH, mostly using a greedy algorithm to iteratively merge nodes, and we adopt the technique of Bradshaw & O'Sullivan (0004). Given the BSH, we propose one for each node pair I and J on two rigid recursive definition of contact potential P puy
7

---

Published as a conference paper at ICLR 2026
bodies, and use the potential of the root nodes final contact potential. Finally, we show that our definition pertains Barrier-Form Smoothness, and Non- Non-prehensile & Non Non-vanishing We start from the base case. For pair of leaf nodes t, U Utj. use Corollary 5.2 to define the following pairwise potential:
d1 RR, + Rí, and d2 = (1 + c) pluu, plus, and P = Plu di P
(6)
di
Specifically, we blend the exact potential between the pair of triangles and the closed-form centered potential, when the distance between triangle centers grows by a factor of e. We leave as user- defined margin that controls the exactness of potential evaluation. Next, given an arbitrary internal node with two child nodes being I and J. we recursively replace the more accurate potential be- tween child nodes with a single potential between parent nodes. Let us suppose I In J Ø, we define potential of similar form as Equation 5
P.n 12 1+
=
which is the potential penalizing distance between two sphere centers. We can only use the centered potential when the two spheres are well-separated, i.e. Ri R, | xy |- Otherwise, we have to use the more accurate potential by descending the tree and sum up the pair-wise terms between each pair of child nodes. Specifically, we define the set of child nodes C(I) [I, if is an internal node and C(I)= (I) I a leaf node. Finally, we define the following potential between internal node:
d1 -R Ry and -2+ pJus= EI,EC(() EIeE((o) P and pJuJ PIuJ
(8)
pdu = di-dd
The main idea behind our formulation is illustrated in Figure 3 and we are ready to present our main result, which shows that the so-defined contact potential well-behaved:
Theorem 5.4. If then E PIUT pertains Barrier-Form Smoothness and Non- prehensile & Non- Von-vanishing) where the summation is taken over the root nodes of different rigid bodies.
This result lays the foundation for our efficient-to-evaluste and well-behaved contact model. Al- though analyzing the cost of evaluating P could be rather difficult for general cases, we follow the idea of fast multiple expansion (Greengard & Rokhlin 1987) and analyze the cost of evaluating the contact potential for a uniform grid, where we show in Appendix th2 the cost is O(T). Finally, notice that contact model should account for frictional contact forces. In Appendix A.3] show that the frictional damping potential proposed in Ye et al. (2025) can be slightly extended to ours.
EVALUATION
We evaluated method in row of five contact-rich manipulation and locomotion tasks: Billiards, Push, Sort, Ant-Push and Gather. We optimize the sequence of con- trol signals using gradient descent at a fixed learning rate minimize user-defined loss functions. More experi- ments and details are deferred to Appendix AA For fair- ness, we compare our contact model with the standard IPC model used in Li et al.) (2020); Huang et al. (2024), and SDRS contact model proposed by Ye et al. (0025). which only violates Non-vanishing We also compare with MuloCo simulator (Tassa et al., 2012). which uses soft contact and provides gradient by finite-difference schemes. Finally, we compare with Suh et al. (2022b) that uses first-order bundled gradient.
Figure 4: Book stacking problem.
8

---

Published as a conference paper at ICLR 2026
Physies Accuracy: Due to method generating contact forces between objects that not in direct contact, we need to validate the physical accuracy of our contact model under different con- tact coefficients μ. We validated the contact model's physical accuracy through the book stack- ing problem (Hall, 2005) in Figure g In this scenario, we stacked 20 planks with dimensions of 2.0m x 0.4m 0.2m and mass of 0.16kg sequentially, extending each plank outward the max- imum theoretical distance from the bottom plank without collapsing. To account for simulation errors, they were shifted inward by 0.1% of the plank's length. We verified whether the system could remain stable under different μ.. Subsequently, we measured the margin between each plank under different μ and calculated the error between the contact force the top plank received from adjacent planks and the theoretical value. The results, shown in Table indicate that the system remains stable when μ < 1e 6, with margin errors in the millimeter range, and the errors between contact forces and theoretical values are negligible.
Contact Coefficient
1e-3
1e-6
1e-7
1e
1e-10 10
Margin (m)
1.82c
5.67e
1.47e -3
4.12c -4
3.13e
Contact force (N)
6.38e
6.40e
6.42e
6.42c /
6.42e /
Success
X
Table 2: The relationship between the margin between the planks, the error in the contact force received by the top plank, and the stability of the system with respect to the contact coefficient p in the book stacking problem.
1000
Dar
600
5345
5045
160
B
#Epoches
3 Epoches
Figure 5: We show the initial (a) and final (b) frame of the billiards task, the convergence history uses trivial initialization (c) and random sampling initialization (d).
Billiards: In this benchmark, have 16 balls the ground with one target red ball whose initial horizontal position and velocity can be controlled. The goal control is for the two green balls to reach the target positions (green circle), where the loss function is the squared distance between the green balls and the center of green circles. We experiment with two different methods for setting initial solutions. Our first method uses trivial initialization where certain rigid objects in scenario are far apart, for which the gradient might vanish except our method. Our second method uses random sampling of control signals to find an initial solution for which gradient information does not vanish. The convergence history of various contact models and initialization strategies is summarized Figure 5 We optimize a trajectory with 100 timesteps at a timestep size 0.04. Except our method, other methods cannot make any progress without sampling due to gradient vanishing, which can be fixed via sampling, while our method achieves faster convergence, with without sampling. We further pick random frames in this scenario and compare the cost of computing the contact potential using our method and IPC in Figure Clearly, our method is slower than IPC due to nested optimization and hierarchical blending. But method accelerated by BSH is much faster than the brute-force version computing all pairwise potentials between triangles.
2013
150
200
Figure 7: We show the initial (a), middle (b), final (c) frame of the push task, and the convergence history (d) our method, IPC-based contact model, SDRS, MuJoCo and gradient bundle method (GB).

---

Published as a conference paper at ICLR 2026
Brute-forca BarteriP
Push: In this benchmark, optimize the po- sition and orientation of a rod to push the red box to reach target red circle, where the loss the squared distance between the box and the red circle. For this benchmark, receding-horizon control to generate a trajec- tory of 200 frames with horizon of only 48 frames, i.e., we iteratively optimize a 488fframe sub-trajes uub-trajectory and only apply the first action. The rod is initially located in front the box, and our contact model guides the rod to first move around to the back of the box, and then push it multiple times to reach the target area as illustrated in Figure 7 This result demonstrates the ability of our contact model to discover multi-stage, contact-rich motions from trivial initial guess. In compar- ison, all other methods make no progress, even with random sampling to find an initial solution that gradient does not vanish. Gather & Sort: In these two benchmarks, we further increase the
Oura
-3.
-6.
et.umberef trinngis
Figure 6: The cost of computing contact potential (left) under different ball mesh resolutions (right).
(6)
Slopa
Figure 8: The final frame of the gather (a) and sort (b) task, with convergence history in (c) and (d).
complexity to have the rod manipulate multiple objects on the table. In the gather task, 10 cubes are disseminated and to be pushed together into the target region. In the sort task, 10 cubes are randomly labeled and to be pushed into separate regions according to their labels. The loss the sum of squared distance between boxes and targets. The frames and convergence history are illus- trated in Figure S Again, other methods make little to no progress, while our method successfully accomplishes the task, relying solely the gradient information. Ant-Push: In this benchmark,
58
100
Figure 9: The initial (a), final (b) frame of the Ant-Push task, and the convergence history (c).
show that our method can work with articulated robot models. We optimize for 9-link, 16-DOF robot to push box to the target position. Once again, our method significantly outperforms other methods, as illustrated int Figure
CONCLUSION
We present detailed analysis of the qualifications for contact model to be well-behaved, which strictly prevents collisions, supports differentiable simulations, induce non-prehensile forces, and avoids vanishing gradients. By hierarchically evaluating the contact potentials assisted by BSH, we further present a well-behaved contact model that is also efficient to evaluate. By analysis on the special case of a uniform grid, we show that the complexity of evaluating our contact potential is linear. Through evaluations on various motion planning and control tasks, we highlight that our optimizer to search for complex motion plans and locomotion gaits that are impossible for previous contact models. Our method is not without its problems. First, we can only handle rigid bodies, and we cannot deal with more general deformable objects for soft robot locomotion soft object manipulation. This is because bound the actual triangles if deformation happens, which could potentially violate Non-prehenssle & Non-vanishing Second, contact potential involves recursive definition and requires nested optimization between pairs triangles, which incurs considerable overhead to a conventional rigid body simulator.
model can guide gradient-
bounding spheres might not
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Pankaj Agarwal, Leonidas Guibas, An Nguyen, Daniel Russel, and Li Zhang. Collision detection for deforming necklaces. Computational Geometry, 28(2-3):137-163, 2004.
Brandon Amos, Ivan Jimenez, Jacob Sacks, Byron Boots, and Zico Kolter. Differentiable mpe for end-to-end planning and control. Advances in neural information processing systems, 31, 2018.
Rika Antonova, Jingyun Yang. Krishna Murthy Jatavallabhula, and Jeannette Bohg. Rethinking optimization with differentiable simulation from a global perspective. In Conference on Robot Learning. pp. 276-286. PMLR, 2023.
Josh Barnes and Piet Hut. hierarchical o (n log n) force-calculation algorithm. nature, 324(6096): 446-449, 1986.
Gareth Bradshaw and Carol O' Adaptive medial-axis approximation for sphere-tree struction. ACM Transactions on Graphics (TOG), 23(1);1-26, 2004.
Tyson Brochu, Essex Edwards, and Robert Bridson. Efficient geometrically exact continuous colli- sion detection. ACM Transactions on Graphics (TOG), 31(4):1-7, 2012.
Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and Zico Kolter. End- to-end differentiable physics for learning and control. Advances in neural information processing systems, 31, 2018.
Asen L Dontchev and R Tyrrell Rockafellar. functiooss and solution mappings, volume 543. Springer, 2009.
Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, and Wojciech Matusik. Diffpd: Differentiable projective dynamics. ACM Trans. Graph., 41(2), November 2021. ISSN 0730-0301.
Susan Fisher and Ming C Lin. Deformed distance fields for simulation of non-penetrating flexible bodies. In Computer Animation and Simulation 2001: Proceedings of the Eurographics Workshop in Manchester, UK, September 2-3, 2001, pp. 99-111. Springer, 2001.
Theodore F Gast, Craig Schroeder, Alexey Stomakhin, Chenfanfu Jiang, and Joseph M Teran. Opti- mization integrator for large time steps. IEEE transactions on visualization and computer graph- ics, 21(10): 103 1115, 2015.
Leslie Greengard and Vladimir Rokhlin. A fast algorithm for particle simulations. Journal of computational physics, 73(2):325-348, 1987.
Eran Guendelman, Robert Bridson, and Ronald Fedkiw. Nonconvex rigid bodies with stacking.
ACM transactions on graphics (TOG), 2713)887 2003.
John F Hall. Fun with stacking blocks. American journal of physics, 73(12):1177-1111, 2005.
David Harmon, Etienne Vouga, Breannan Smith, Rasmus Tamstorf, and Eitan Grinspun. Asyn- chronous contact mechanies. ACM Trans. Graph.. 28(3), July 2009, ISSN 0730-0301. doi: 10.1145/1531326.1
htl.http:: /tttss://os. doi.. org/10...0.. 01445555132665 1
1531393
Eric Heiden, David Millard, Erwin Coumans, Yizhou Sheng, and Gaurav S Sukhatme. Neuralsim: Augmenting differentiable simulators with neural networks. In 2021 IEEE International Confer- ence on Robotics and Automation (ICRA), pp. 9474-9481. IEEE, 2021.
Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua Tenenbaum. William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik. Chainqueen: A real-time differentiable physical simu- lator for soft robotics. 2019 International conference on robotics and outomation (ICRA), pp- 0265-62711 IEFF, 2019.
Zizhou Huang, Davi Colli Tozoni, Arvi Gjoka, Zachary Ferguson, Teseo Schneider, Daniele Panozzo, and Denis Zorin. Differentiable solver for time-dependent deformation problems with contact. ACM Transactions Graphics, 43(3):1-30, 2024.

---

Published as a conference paper at ICLR 2026
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6990, 2014.
Simon Le Cleac'h, h, Mac Schwager, Zachary Manchester, Vikas Sindhwani, Pete Florence, and Sumeet Singh. Single-level differentiable contact simulation. IEEE Robotics and Automation Letters, 8(7):4012-4019, 2023, doi: 10.1 109/LRR..022333688822
Quentin Le Lidec, Igor Kalevatykh, Ivan Laptev, Cordelia Schmid, and Justin Carpentier. Differen- tiable simulation for physical system identification, IEEE Robotics and Automation Letters, 6(2):
1413-3420, 2021.
Sergey Levine and Vladlen Koltun. Guided policy search. In International conference on machine learning, pp. 1-9. PMLR, 2013.
Minchen Li, 7achary Ferguson, Teseo Schneider, Timothy R Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M Kaufman. Incremental potential contact: intersection- and inversion-free, large-deformation dynamics. ACM Trans. Graph., 39(4):49, 2020.
Sizhe Li, Zhiao Huang, Tao Du, Hao Su, Joshua Tenenbaum, and Chuang Gan. Contact Points Discovery for Soft- Soft-Body Manipulations with Differentiable Physics. In International Conference on Learning Representations (ICLR), 2022a.
Yifei Li, Tao Du, Kui Wu, Jie Xu, and Wojciech Matusik. Diffcloth: Differentiable cloth simulation with dry frictional contact. ACM Trans. Gruph., 42(1), October 2022b. ISSN 0730-0301.
Zhehao Li, Qingyu Xu, Xiaohan Ye, Bo Ren, and Ligang Liu. Difffr: Differentiable sph-based fluid-rigid coupling for rigid body control. ACM Transactions Graphics (TOG), 42(6):1-17, 2023.
Chen Liang, Xifeng Gao, Kui Wu, and Zherong Pan. Second-order convergent collision-constrained optimization- planner. IEEE Robotics and Automation Letters, 2024.
Pingchuan Ma, Tao Du, Joshua Tenenbaum, Wojeiech Matusik, and Chuang Gan. RISP: Rendering-invariant state predictor with differentiable simulation and rendering for cross-domain parameter estimation. In International Conference on Learning Representations, 2022. URL
https: Btttp::/openreview net wttp://wwwwweettotteeeeeeme
Jerrold E Marsden and Matthew West. Discrete mechanics and variational integrators. Acta numer- ica, 10:357-514, 2001.
Igor Mordatch, Emanuel Todorov, and Zoran Popović. Discovery of complex behaviors through contact-invariant optimization. ACM Transactions on Graphics (ToG). 31(4):1-8, 2012.
Rhys Newbury, Jack Collins, Kerry He, Jiahe Pan, Ingmar Posner, David Howard, and Akansel Cosgun. A review of differentiable simulators. IEEE Access, 2024.
Zherong Pan and Dinesh Manocha. Active animations of reduced deformable models with environ- ment interactions. ACM Transactions on Graphics (TOG), 37(3):1-17, 2018
Tao Pang, HJ Terry Suh, Lujie Yang. and Russ Tedrake. Global planning for contact-rich manip- ulation via local smoothing of quasi-dynamic contact models. IEEE Transactions robotics, 2023.
Clemens Schwarke, Victor Klemm, Joshua Bagajo, Jean Pierre Sleiman, Ignat Georgiev. Je- sus Tordesillas Torres, and Marco Hutter. Learning deployable locomotion control via differ- entiable simulation. In 9th Annual Conference on Robot Learning, 2025.
Tuur Stuyck and Hsiao-yu Chen. Diffxphd: Differentiable position-based simulation of compliant constraint dynamics. Proceedings of the ACM on Computer Graphics and Interactive Techniques, 6(3):1-14, 2023.
Hyung Ju Suh, Max Simchowitz, Kaiqing Zhang, and Russ Tedrake. Do differentiable simulators give better policy gradients? In International Conference on Machine Learning, pp. 20668- 20696. PMLR, 2022a.
12

---

Published as a conference paper at ICLR 2026
Hyung Ju Terry Suh, Tao Pang, and Russ Tedrake. Bundled gradients through contact via random- ized smoothing. IEEE Robotics and Automation Letters, (2):440004000) 2022b.
Yuval Tassa, Tom Erez, and Emanuel Todorov. Synthesis and stabilization of complex behaviors through online trajectory optimization. In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 4906-4913, 2012.
doi: 10.1109/IROS 20w.H/..........86..// 2.686620.
Emanuel Todorov. A convex, smooth and invertible contact model for trajectory optimization. In 2011 IEEE International Conference on Robotics and Automation, pp. 1071-1076, 2011. doi: 10.1109/ICR..011. 1.5999144.
Mare Toussaint, Kelsey R. Allen, Kevin A. Smith, and Joshua B. Tenenbaum. Differentiable physics and stable modes for tool-use and manipulation planning extended abtract. In Proceedings of the Eighty-Eighth International Joint Conference Artificial Intelligence, IJCAI-19, pp. 6231 6235. International Joint Conferences on Artificial Intelligence Organization, 2019. doi: 10. 24963/ijcai.20 201..200996666 httpsss:://aai wrg/10
or//:///////....0/..././. 22000249966 249633314aa1.211177855 2019//66
Dylan Turpin, Liquan Wang, Eric Heiden, Yun-Chun Chen, Miles Macklin, Stavros Tsogkas, Sven Dickinson, and Animesh Garg. Grasp'dt d: Differentiable contact-rich grasp synthesis for multi- fingered hands. In European Conference on Computer Vision, pp. 201-221. Springer, 2022.
Keenon Werling, Dalton Omens, Jeongseok Lec, Ioannis Exarchos, and C. Karen Liu. Fast and Feature-Comples Differentiable Physics Engine for Articulated Rigid Bodies with Con- tact Constraints. In Proceedings of Robotics: Science and Systems, Virtual, July 2021. doi:
Jie Xu, Tao Chen, Lara Zlokapa, Michael Foshey, Wojeiech Matusik, Shinjiro Sueda, and Pulkit Agrawal. An End-to-End Differentiable Framework for Contact-Aware Robot Design. In Pro- ceedings of Robotics: Science and Systems, Virtual, July 2021. doi: 20.15607/RSS.222.XXV 008.
Jie Xu, Miles Macklin, Viktor Makoviychuk, Yashraj Narang, Animesh Garg, Fabio Ramos, and Wojciech Matusik. Accelerated policy learning with parallel differentiable simulation. In In- ternational Conference Learning Representations, 2022. URL http:sss//ppeeeeeeleee
nett/Eorum:1d==Z3KR0M
Xiaohan Ye, Xifeng Gao, Kui Wu, Zherong Pan, and Taku Komura. Sdrs: Shape-differentiable robot simulator. IEEE Transactions Robotics, pp. 1-20, 2025. doi: 10.1109/TRO.202536666644
13

---

Published as a conference paper at ICLR 2026
A APPENDIX
A.1 ADDITIONAL PROOFS
proof of Lemma 4.1
Barrier-Form If æ Crose then each pair of t, t, on different rigid bodies is disjoint as the convex hulls are closed sets. Therefore, by the separating plane theorem, there exists a separating plane Psy and some positive Ciz > 0 such that:
|mee|l 1 (II(())EL) m4y)) di e2j/2 (I)())tt, mij) dij S -eig/2.
Clearly, Plotts (2s(P::/2, { co. On the other hand, if Cabs, then there exists non-disjoint pair t, for any separating plane Pij there exists (Z.(k): (Z)(k)))B))) dij 0 or Ty(1), naj) + dig O, leading plut, co. Finally, at any feasible solution Pij. we must have [|ny > because otherwise, have:
Cij
Σ
co.
(dg))
(-diy) k=1
We have thus established Barrier-Form for each plut and thus P.
Smoothness This follows from the inverse function theorem (Dontchev Rockafellar 2009), the smoothness of problem data Lij. and strictly convexity of Lij.
Non-prehensile Non-vanishing We can simply define A(x) = ( (ti,tt))4 # {2} such that every pair of disjoint triangles ty) appears in exactly one term of plut, Thus, we only need to verify that stattts ((kk)t, e Ft,,,,, and the case with j(K)et, is symmetric. By the implicit function theorem, we can derive the analytic formula:
douts
Aij ((F:()) + dzy)²' i(k)et,
CH())(H)e))) and & > 0, then we have miy) + d Fey-t,- Therefore, we can in turn prove the sufficient condition that n. = a(a - b). Due to the optimality of Lis with respect to Pij. we have:
Suppose ni = a(a - b) for some a E CH(II(N))t,)) (((ket,; a(a
OLij
nij
Fi(k)
T((h)) ((15(K) mii)) day)
12
Onij |ni
(ITA)) )
R-
0hg
3
Σ
Odij
(ITE()) miy) )z ((I)(K) niis) d
From the above two equations, we can conclude that nij aaa b) by defining:
||ncs| (1
||nss
3
1
k=1 (Ia(N),ne)) dig)² 00
12
3
Ti(k)
Σ
Re (ZZ(R) R(y) ) di3 e CH((,(())))
a= =
I (I( miy) di3) )2
kz(k))
3
1
e &CH(E)()))))
=
1= (1)(()) -nip) ) diz)²
k= ((25(K), -nis) d
thus all is proved.
Proof of Lemma 5.1 Barrier-Form Case I: When || - < di, PTuj = which per- tains Barrier-Form since p4u. pertains Barrier-Form by assumption. Case II: When || - x5| d1 2 Rz R,,, there are two sub-cases. Case IL.:: If r Cos e.g. where two bounding spheres just touching and the touching point lies common triangle, then 0 PIuJ PTuJ where the first two inequalities are due to our assumption and the last equality due to PluJ
pertaining [Brriir-Form] PIUJ PTuJ 00. Case II.b: If CIUJ then have
di-d2
14

---

Published as a conference paper at ICLR 2026
0 PTus PT.. < following the same reasoning case ILa, 0 S P7 < co. Case
III: When |azz - Σzll
d1 2 RI + Ry, then we must have x CTUJ and the same analysis as case
co. Thus, we have verified that P pertains Barrier-Form in all cases. I1.b leads to 0 di ddg
Smoothness This due to Smoothness in P PTuj, and the second differentiability of Qd,-da-
Non-prehensile & Non-vanishing The force Fez induced by di takes the following form:
8phus
IuJ = -(1 - Qd -dy (a)) Uzut
Od -dg (x)
0zie
torm
torm 11
d'(112x 2J d1) / di))
(PPt -Pau
(II - (z)-
(d2 d1 21 - x 111
There are three terms and we show that each term belongs to FI-II. For term I, we know that -OP2.0 /0xzz E FI-I since d, satisfies Non-prehensie & Non-vanishing Since the coef- ficient (1 odd (x)d((x) 2 0 and Fイ-Z is cone, we conclude that term is zero belongs to FG-I. The same reasoning applies to term II. For term III, 3I - Iy belongs to FJ-I. Since plut s ppIuJ by our assumption, the remaining coefficient is non-negative, thus term III is zero do
belongs FI-I. Finally, least one of term I or term II is non-zero, so we conclude that fluJ eFy-I and all is proved.
We first show that Equation 5 is correct. The objective function in Equa- tion is derived by replacing all Tx(k) and T)(K) in Lij with the center points Ft, and Fty, respec- tively. By symmetry, the optimal separating plane must be the middle surface between Dt, and Zty>
Proof of Corollary 5.2
taking the following form:
nsy a - II,)) and di = -a/2(zz,, -I, xe,)
Plugging the optimal separating plane and solving for leads to Equation | Next, we show that all three assumptions in Lemma [.1 hold. In fact, there are only two non-trivial assumptions. We first
show that plu pr.tt when | æ d1. Let use denote piy as the optimal separating plane for pl.ut; then we have the following inequality:
1
Plule argmin
12 (1
(nil))* I)
(TTe, + dig)*
k=1 (1Te, -mi))
) d
Ps
k=
1
<12 (1
l)
((Z,, miy) + +dz))=
(1Ze,) -miy)) dii))*
k=1
<12
(1 || l))
k-1 (1i2a) m,y) + dis)
hm ((zz()) -mty) - diz))
-ptuts
where the first inequality is due to optimality of P.ut the second inequality is due to the convexity of function 1/((e, mij) d and 1/((o,-n)) -n;;) d dr We then show that phut, satisfies Equa-
tion 2 The force any T4(k) takes the following form:
4(xr, Fz,))
1
JUJ
It. Iz, 15/2
|0,, It, 1/
which clearly belongs Fz thus all is proved.
Lemma A.1. If'e 0 then PIuJ pertains
Barrier-Formlfor any node pair I, of two rigid bodies.
Proof. First, by induction from leaf to the root node, we can verify that 0 s PIuJ BSH Second, suppose
E Crree, then all the wiie terms between leaf nodes pfutt; Further, for all the centered
15

---

Published as a conference paper at ICLR 2026
potential in Equation 7 we have Plus < because they evaluated only when Ry Rg | Iz I|. The root potential PASH is then derived by finite number of blending and summation so we have P co. Third, any x Cobs. there exists non-disjoint pair 4, Lj belonging the two rigid bodies. We will show the following two claims hold by induction from leaf to root:
oo at any node such that 4 S I and 1, S.
BSH
Base Step: BSH
Pai = oo by Corollary 5.2
Inductive Step: We assume our first claim holds for any t, E Ie E C(I) and t, Jo C(J). If and tj S JJ we must have laz - -I5ll Rt Ry due the pair 4,tj being non-disjoint. The children set satisfying & E Ie e C(I) and t, é Je é C(T) can always be found, so we have plu. =PZ IpTuu 00.
BSH
To prove that P pertains Non-prehensile Non-vanishing i.c., the non- non-prehensile and non- vanishing property, also need to induction. To this end, establish the non-prehensile property for an index subset as follows:
Definition A.2. The pairwise potential PISH pertains Non-prehensile Non-vanishing restricted to (I,J) if. at every I E Cree. we can define finite family of set pairs AzuJ (x) such that RSH E(I' J'eAAs PESH where every term pI''u3' satisfy Equation Further, for every pair of disjoint triangles (t,t,) such that t; and t5 I or vice versa, we have t, Lt, = least one (I' .J') 6 AAuu (z).
for at
Lemma A.3. c O then BSH pertains Non-prehensile & Non-vanishing Nonn-vanishinggrestiicce restricted to (I,T) for any node pair I.J of two rigid bodies.
Proof. At x Crree, we show the following two claims by induction from leaf to root:
The pairwise potential PZUJ
P.UJ
d2
at any node when || xz Rr Ry.
The pairwise potential (I, T) for any node pair.
BSH
pertains
Non-prehensile & Non-vanishing restricted to
Base Step: For the pairwise potential Plutt = -dd we define Acque, (x) ((ti, {, then Non- prehensile & Non-vanishingaan and the fact that BSH 2 p!lut follows from Corollary 5.2
Inductive Step I: We assume first claim holds for all I. C(I) and Je E C(J), i.. DI.UTe
PZ..U. when [zz xTJ_l Rr,- R5g. We first show that PIuJ pIuJ when [ xJ Rr Rg. We note that I J I contains at least 3 triangles, otherwise we reduce to the base case, so there least 2 terms of form PI.UJ. Each such term has the following lower bound:
P.uJ. 2 Plinno 12
1
2 12
1+
II.
Igg
Rx + R
Here, the first inequality is due to inductive condition and the fact that BSH is a layered hierarchy by Definition A.2 so that | x 2 Rz + Ry implies ITz Izz 2 R.g. + R. The second inequality is because XI, (resp. x5.) is most Piz (resp. R2)) from Fr (resp. xg). )) Using
16

---

Published as a conference paper at ICLR 2026
the above lower bound, we derive the following estimate:
2
Paus 224 1+
+ Rr + RJy
48
24
224 +
v |
IJ
+ RR Rg
RI + R3
48
24
224 +
+
/2I PI 35
2 1z xz
2
212 +
24
12
12
1+
1
plus = -plu
v || - xzl
| Izl
Fzl
where we use the fact that 1|oz - J.5l I 2 d1 Rx + Rg in the third inequality. As a result, we have first claim holds for Iu J by the definition of the blending Equation
Inductive Step II: We assume second claim holds for all I C(I) and Je e i(e), P I.u.Ja pertains Non-prchensile & Non-vanishinglrestrieee
nsu
Jo). We show that PLyI pertains Non- prehensile & Non- vanishing restricted to (I.J) by considering three cases. Case II.:: If 4 ا S d1, then Pasi Paus consists of terms of form Pisu each satisfying Non-prehensile & Non-vanishing by our inductive condition. Let us now define the finite family of set parrs by (Ir). It can be verified that this union is a disjoint union, and for every pair of disjoint triangles (1: E I, Lg J), we have 4, Lu I'' u J' belongs to exactly one such AI,uJ, (z)). Further, we have PITSH
the union Azur (z)
PT"'u", where each PTu'" satisfies Equation 2 due to Non- & Non-vanishing of the corresponding PESH P...J= by our inductive condition. We have thus verified Non-prehensile & Non-vanishing of BSH Case II.b: If || xz 2 d2, then P = -pZu p.ud is singled, centered potential. We trivially define Azur (x) [(I then clearly every pair of disjoint triangles t, Ut, IIUJ (Azu(()). Further, the induced force takes the following form:
Zu.J
12(xg - xz)
IUJ NeT
1+
||| x 3//2
|
xg
1/2
which clearly belongs to F.-II, thus we have verified all conditions in Non--rehensile & Non- vanishing of P I.JJ Case II.e: If < <|>z-zz|| < d2, then PESH is blending of P in Case IL. and Ploy Case II.b, with strictly positive weights. We again trivially define Azur (z) = ((I.T) T)) to have every pair of disjoint triangles 1, L S IuJ e AAAu(()) The only condition we need to verify is that PTTSH satisfies Equation 2 as single term. We use a similar technique as in the proof of Lemma by expanding the force term to get Equation @heeee there are three terms. For term II, We have by the analysis in case II.b that - /@xiez 6 FFJI and Φd (I) > 0 is strictly positive, so term belongs to Fy-zz Term III is zero or belongs to F3-I because DIUJ PTuJ by our first claim. For term I, we know from the analysis Case II.a that:
BSH
JpZus
IpI'
Σ
BSU
Ox
dx,et
(I' (I)
where each term -OpT'uJ' BSH 10xHI is zero or belongs FJ F Combined with the fact that the coefficient (1 - Pdi (x)) > 0 is strictly positive, we conclude that term I is zero or belongs
to F As a result, we see that Tuu Non ranishing
satisfies Equation as single term, so Non-prehensile &
Non-prehensile & Non- vanishing follows from Lemma A. and Lemma A.33 respectively. Smoothness follows from the fact that P is derived by finite number of blending between pairwise potentials, and all potentials and blending operators are twice differentiable.
Proof of Theorem 5.4 Barrier- Form and
17

---

Published as a conference paper at ICLR 2026
Level 2
Level
Level
Level-1
Figure 10: A special BSH constructed for the uniform grid. The levels of the BSH are indexed from bottom up. The 11h level contains only leaves. Every (2i + 1)th level (i 0) merges two diagonal, rectangular blocks. Therefore, every (2i)th level consists of a cubic mesh block of side length 2
A.2 COMPLEXITY ANALYSIS FOR UNIFORM GRID
We show that, in the special case of two rigid bodies in the shape a square uniform grid with infinitesimal distace each other, as illustrated in Figure 11 the cost of evaluating PESH is O(T). Here we assume a 2D uniform grid with N2 grid cells so that generality, assume the grid size of BSH as illustrated in Figure TO It is easy to see that the bounding sphere radius of nodes in each level all the same, which is denoted Tt- We have the following results for estimating Ti which can be derived directly from the construction of PZUJ.
Further, for ease of analysis, O(N3). Without loss of adopt special construction
T-1 V5/3 TO 22+22/5
72:-1 = $2 1/2, Vi O with C, I+vI0
To derive the second property, note that the tightest bounding sphere for the geometry of a (2i)-level node is 2/-1/2 This implies that the tightest bounding sphere for the geometry of 0-level node is 1/v2. However, the actual = (VT+ /10)/3, so the bounding sphere is unnecessarily scaled by C,. By induction, we can prove that we can scale all tightest bounding spheres by C, accordingly to satisfy Definition 5.3
TuJ To analyze the cost of evaluating BSH we note that, the cost reduces evaluating series of interaction terms either using centered po- tential (PIuT) or between leaf nodes (Pt.tt) Further by the recursive definition and the spe- cial structure our BSH, the interaction terms always computed between two nodes at the same level. Therefore, we can upper bound the number of interaction terms level by level.
T
Case I: (2i)-level Node We focus on the nodes at (2i)th level, which bounds a cubic mesh block of side length 2 We can de- note these mesh blocks as Bi indexed us- ing subscript
Figure 11: Two rigid bodies in the shape of square uniform grid with infinitesimal distance d d O.
In other words, B sists of all triangles with coordinates within [m22, (m 1)2') [n22, (n 1)2']. For two
18

---

Published as a conference paper at ICLR 2026
such blocks indexed by B I and B I on the two rigid bodies, can define their distance as:
(BB B n max(Im m n'l).
Now suppose d(B man B m'n'') 2 2, we can be sure that E and B belong to different (2i 2)-
level blocks. Specifically, we define B s relationship:
B2i=2 mn and [B
s [2:+2 and we have the following m/w
d(BBa
B mmi
d( [ 2i+ mn
(10)
2
Now let us define I and J mesh blocks is sufficiently faraway, the potential term pzud reduces to the centered potential Pluz be computed via Equation / without utilizing any information from lower levels. A sufficient
2i+2
m/Rt we know that the distance between the two cubic
condition for this to happen is as follows:
d( 21
B )2' 2r2i+2(1 +e.
(11)
Combining Equation 10 and Equation 11 we know that if d(B [ Bmx 2
then the interaction between the two (2i)-level blocks would be handled by the two (2i level super-blocks. Therefore, in order to compute the contact potential, we only need to evaluate the interaction between 2 and at most (1 4/2v/2, e)2 blocks around it. The number of (2i)-level interaction terms of form pTuJ is:
ol N (1 (1+412V5,,(++0)))
where the first part [ is the number of (2i)-level level blocks and the second part is the number of other (ei)-level blocks, with which an interaction term pluJ needs to be calculated.
Case II: (2i - 1) leeell Node We have finished analyzing the cost of (2i)-level node interactions. The case with (2i level node iterations (i 1) almost identical. Each (2i 1) level node consists of half the triangles of (2i)-level node, so we can assume the (2i)-Ievel node B has two children denoted as Fnn and BBi--... Without loss of generality, we only consider mn which has the same bounding sphere center as that of B By the analysis of Case we know that, for all the children of B with d(B22 (BB [B2 ( 2 2[2v/C, + <)], their interaction with B2i-1d would be taken care at level (2i + 2). Again, we conclude that we only need to evaluate the interaction between BB2--. and the children of most (1 + 422 (1+ee,,,++1)))))) (1+ (2i)-level nodes
2
me around The number of (2i 1)-level interaction terms of form DTuJ is again:
o ((
N 2
(1 (1+412/2C((L+))))
Case III: 1-level Leaf Node The case with leaf nodes exactly the same as that (2i- 12i-ve-level nodes. Each O-level node has two children -1-leeel denoted as B and B Focusing on B and by the analysis of Case II, we know that, for all the children of B with d( (BB) B",,.)) 2 2 11,22,,(1 (1+ e. their interaction with B would be taken care of at level 2. Therefore, conclude that we only need to evaluate the interaction between B- and the children of at most
(1+ 4[2/2C, + +) O-level nodes around it. The number of (2i level interaction terms of form pl.ut is again:
Put everything together, the cost of evaluating PASH is:
[oga N1
N
0(N3) = =
200
19

---

Published as a conference paper at ICLR 2026
A.3 FRICTIONAL CONTACT MODELING
A feature-complete potential function should be able to handle frictional contacts. To this end, we adopt the technique proposed by Li Lt al. (2020) Ye al. (2025) and consider the contact potential between the pair of triangles ( 2 (2) (3) (1)) & (()) (()))>((3))
(33))) where we write explicitly the six vertices related to this potential function. The negative gradient norm of this potential is the normal force applied on the corresponding vertices. [Li al.] (2020) proposes to model the contact
potential as tangential velocity damping term weight by the normal force magnitude. Put together, the frictional damping term takes the following form:
3 Σ 1=1
plut, (i(k)
pluee ((K)
D(x 11 x., 5t)
Σ 1
Pi((((()) (+1 Dii
((K) ((6) 5t) + Σ A
D
where the summation is taken over triangle pairs on different rigid bodies, Dj is the tangential velocity damping term, which penalizes the relative velocity between 4; and 1j. and À is the frictional coefficient. However, the above formulation is not strictly second-order differentiable. To fix this problem. Ye al. (2025) proposed novel definition of Dj. which does not penalize the relative velocity between it and Ej. Instead, they penalize the relative velocity between the two triangles and the separating plane Pij. assuming the plane a physical object with zero-mass. We refer readers to their work for more details.
The original frictional damping term can only deals with frictions between triangles. However, our contact potential DT. is a hierarchical blending of potentials between triangles and centered potentials. To extend the frictional damping term to use our PTuJ centered potentials and only consider potentials between triangles. Specifically, we propose the following potential between a pair of triangles:
BSH we propose to disregard the
d1 R2 + +R, and d2 = (1+e)d1 P4i plue, and pluut pliut du-da
(12)
local
Compared with the potential between leaf nodes in Equation Equation 12 is not globally supported and vanishes when the distance between triangles larger than d2, denoted using subscript local- This design choice would not cause gradient vanish because our normal contact potential P always provides non-vanishing gradient. In parallel, the benefit using a locally supported function that we can use a bounding volume hierarchy to quickly reject faraway triangles as done in [i et al. 20200;; Ye et al. (2025). Equation formulation:
12 is plugged into the frictional damping term to yield our final
P fut, local x ((())
plaut ((()))
D(*tt1 t+1 (t,8))== =
Σ b
Σ 1=1
Dj Di(att ((K) at (t)))t)) 0t) + Σ A k=1
Di (h
Intuitively, our frictional damping model assumes that two triangles can only impose frictional damping forces each other if their distance is less than d2- Otherwise, the two triangles can only impose normal forces, but not frictional damping forces. This design choice preserves the property of non-vanishing gradient, but also allows efficient evaluation. Finally, emphasize that as μ 0, both our contact potential and frictional damping term converges to the exact frictional contact model with Coulomb friction.
20

---

Published as a conference paper at ICLR 2026
A.4 SPICE-DIFERENTIABIIE OF IPC
We show that the IPC contact model (Li ct al. 2020) is differentiable but not twice differen- tiable. Let us consider simple 2D case, where the only geometric primitive pairs that incur collision potential is between a point and a line- segment. Let with a single geometric primitive: a line seg- ment with two end points located (1,0) and (2, 0) and point moving on the line (x, 1) with x e [0,3]. The IPC potential for this toy example is formulated as: P(x) - ll))(()))) where d(x) is the distance between the point and the line-segment, parameterized by sin- gle scalar æ. Clearly, the differentiability of P relies on the differentiabilityy of d(x). As ana- lyzed Li et al. (2020), d(x) is differentiable, Od( /0x is well-defined. In Figure 12 plot the value of ad(())/Or when [0,3]. Clearly, Bd(x)/Ox is not a differentiable function, so we conclude that P(x) cannot be twice-differentiable
0.5
0.0
assume the toy example
-0.5
2
3
X
Figure 12: We plot the value of Od(x)/ /dx when e [0.3] in toy example under the IPC contact model.
The non-smoothness is due to switching between different Voronoi regions. When 1 2 2, the closest point on the line-segment is a vertex. Instead, when x E (1, the closest point lies interior to the line segment.
A.5 HIERARCHICAL POTENTIAL BLENDING
As the same in al. 2020 Ye al. (025 use the Newton method equipped with line-search to solve the optimization problem of Equation where we recursively accumulate the contact potential, its gradient, and Hessian based on BSH mentioned in Section to accelerate the contact potential evaluation. Specifically, we establish a BSH for each joint in the scene and then recursively perform collision detection based on BSH while accumulating potential. When detecting potential collision pair where nodes I and J are both leaf nodes, the distance between two nodes is greater than d2, it reaches the end of the recursion and return PJuJ and d2 respectively. Otherwise, we continue to search for their child nodes and blend it's potential using Equation We summarize the algorithm of hierarchical potential blending in Algorithm
Algorithm Hierarchical Potential Blending
function PROCESS PAIR(I,T)
pluJ +
p + -IIUJ > Equation 7
if Tz xz
d2 then
BSS Pdz -II.
return
if I is leaf and J ty is leaf then PIuJ pliuts D Section
Plug Pdj Pd2 in Equation for Pdj-1d2
PIUJ -PJ
return
for each pair In, J. in < IIJ >'s children do
PROCESS PAINSS__AAR(e..
DIuJ + = P...J.
BSH
Plug d, Pd, Equation for P
pIuJ pIuj
BSSH
di
return
21

---

Published as a conference paper at ICLR 2026
A.6 EXPERIMENTAL DETAILS
In this section, we provide experimental details and extended evaluations.
Parameter
Billiards
Push
Ant-Puse ant-Push
Sort
Gather
Gather Gather-Bunny
trajectory horizon H
100
200
240
700
300
550
receding horizon h
48
16
32
32
potential coefficient μ
1e-7
1e-
le-5
3e
3e-6
5e-
& learning rate Adam (31, B2)
3e-2
le-2
3e
1ea
le
le
(0.3,0.5)
(0.3,0.5)
(0.3,0.5)
(0.3.0.5)
(0.3,0.5)
(0.3,0.5)
number of iterations
400
50
100
100
60
60
degrees of freedom timesteps At
96
12
22
66
66
66
0.04
0.04
0.02
0.04
0.04
0.04
Parameter settings for different benchmarks. Here, the trajectory horizon H represents the total number of frames in the entire trajectory. For benchmarks requiring receding horizon execution, the receding horizon h represents the number of frames in the sub-trajectory. The potential coefficient μ represents the contact energy coefficient, where smaller y indicates more physically accurate contact mechanism. a, the learning rate, denotes the step size for optimization, and the number of iterations specifies the total optimization steps. For all the benchmarks, we set the and 88, to small values, which helps escaping from local minima. At represents the timestep for each frame. Our simulator allows stable, penetration-free simulation even under relatively large timesteps.
Table 3:
of the Adam optimizer (Kingma & Ba 2014), $1
A.6.1 BASELINES
We compare our method against three state-of-the-art baselines. The first the IPC contact model (Li et al. 2020). employed in the differentiable simulator of (Huang al. 2024). How- ever, due to the model's lack of twice-diffcrentiabiiiy the implicit function theorem does not apply. As a result, Huang al. (2024) resort to a few iterations of gradient descent to approximately solve Equation| While this yields usable gradient information, it is well known that the gradients can vanish when the interacting geometric primitives are not in close proximity. The second baseline is SDRS contact model (Ye et al. (205), it's a differentiable model with twice-differentiabbliy But similar to IPC, the gradient also vanish when the interacting convex hulls are far apart. The third baseline is the Gradient Bundle (GB) method (Suh et al.) 2022b), which addresses gradient vanish- ing through sampling, evaluated in practice via Monte Carlo methods. However, when primitives are far apart, the likelihood of sampling non-zero gradient decreases significantly. Consequently, gradients obtained from GB can be both noisy and prone to vanishing with high probability.
A.6.2 BENCHMARK DETAILS
We implement a full-featured rigid body simulator based on our novel contact model. Each bench- mark scenario includes controlled objects and target objects. The controlled objects are actuated using a built-in and contacts to move the target objects to their designated spatial positions. In all experiments, we use the following loss to measure the progress of optimizers:
controller, and the objective across all benchmarks is manipulate collisions
ReLU([ - =comll -
where FCOM and #COM the position and desired position of the center of mass of some target object. étargot is the error coefficient, indicating that tasks are considered successful for certain objects when they are within Enaget of the goal. The statistics of benchmarks are summarized in Table 3
Billiards In this task, the indices of the two target balls and their target locations are randomly se- lected. The objective is to control a distinct red ball that it strikes the target balls through contact, moving them to their respective target positions. Each ball has 6 degrees of freedom, resulting in total system of 96 degrees of freedom. We only control the initial horizontal positions and velocities of the red ball, corresponding 4 control dimensions.
Push In this benchmark, the task goal is to control a rod to push a box to the target region. The system consists of 2 objects with total of 12 degrees of freedom. At each timestep, continuous
22

---

Published as a conference paper at ICLR 2026
6-dimensional control signal is generated to control the rod. The control signal each timestep is obtained by solving a receding-horie -hoiinnn optimization problem.
Ant-Push In this benchmark, our goal is to drive the ant robot to move and push the box to the target position. The ant consists of a base, four large legs, and four small legs. The base includes 3 translational degrees of freedom and 1 rotational degree of freedom, while the upper legs and lower legs connected using ball joints and hinge joints, respectively. As result, the ant has total of 16 degrees of freedom in the kinematic state. Combined with the 6 degrees of freedom of the box, the system has a total of 22 degrees of freedom in the kinematic state, of which we can control 12 degrees of freedom in the ant's legs. We use accumulated sine wave signals to parameterize our controller for each degree of freedom of the ant's legs done in Hu al. (2019). In this case, the decision variables are the amplitude, frequency, and phase of the sine waves.
Sort In this benchmark, two types of cubes are mixed together on the ground. Use a rod to push each type of cube to its target location without mixing with each other. We set up 10 target cubes and one rod, so the system kinematic state has a total of 66 degrees of freedom and we can control the 3 translational degrees of freedom of the rod via PD controller.
Gather In this task, 10 cubes are randomly put on the ground. Use a rod to collect all the objects into one area. Again the task has a total of 66 degrees of freedom and we control 3 translational and I horizontal rotational degrees of freedom of the rod. To further validate the efficiency of our method, we can handle objects with more complex geometries in this scenario. We replace the cubes with bunnies and successfully complete the gather task, we call this benchmark Gather- -unney
A.6.3 ABLATION STUDY OF CONTACT PROPERTY
In Figure 3 we discuss the various properties that a contact model needs to possess. Here, we will analyze the impact of these different properties on the contact model's performance in simulation and policy learning. Collision models without Barrier-Form would allow intersection between ob- jects. It is well known that the distance function between objects is non-smooth when the objects are in collision. Therefore, if we take away |Barrier-Fore Form then [Smoothness] will be violated auto- matically. Non-prehennile requires that normal collision forces between objects are pushing them apart, instead of pulling them together. Almost all existing collision models satisfy this property. We conduct an ablation study on Non-vanishing Since the only difference between our method and SDRS model (Ye et al. (0255) is the additional satisfaction of Non-prehensile] the results of this ablation study can be observed in the comparison between our method and SDRS across vari- ous benchmarks in Section 6 This demonstrates the necessity of satisfying this property for policy convergence when objects are far apart.
We note that it is possible to take away Smoothness alone and conduct an ablation study. To this end, we modify our Equation 3 to use Li et al. (2020) instead of Ye et al.) (2025). In other words, we combine Li et al. (2020) with a globally supported barrier function instead of the original locally supported version. By doing so, we still have Non- vanishinnihggggee but fails [Smoothness] We tested the results the Gather and Push benchmarks in Figure [3] showing that without [Smoothness) the policy convergence speed is significantly slower or may even fail to converge completely.
Ye
A.6.4 INFLUENCE OF STIFFNESS AND TIMESTEP
It is well known that the stiffness of a system and the timestep Đl value can affect the evaluation of gradients. Here, we validate the sensitivity of our contact model to stiffness and timestep. We demonstrate the impact of stiffness and timestep on policy optimization by verifying the conver- gence of strategies on the Push task. our contact model, stiffness is determined by the contact coefficient μ4, where a smaller j corresponds to a stronger system stiffness. Firstly, we fix l 1 and record the simulation time required for strategy convergence under different timesteps. We ob- serve that strategies converge faster with larger timesteps. This is because, when using MPC for decision-making a larger timestep results in fewer frames for the same horizon length, reducing the number of gradient multiplications and avoiding the issues of gradient explosion and vanishing gradients that are common in optimization problems. Since our method is unconditionally robust to timestep size, for a specific task, we should choose the largest timestep within reasonable range to
23

---

Published as a conference paper at ICLR 2026
+
W/O Property 3.2
4
Property
3
4
Logs
Loss
2
2
1
W/o Property 3.2
0
Property 3.2
50
100
150
200
100
200
300
# Steps
# Steps
Figure 13: The convergency history with or without Smoothness Push (left) and Gather (right) benchmark.
accelerate strategy convergence. Secondly, we fix öt 0.04s and record the simulation time required for strategy convergence under different values. We find that only when the system stiffness is ex- tremely high (μ4 1e-19), leading to drastic gradient changes caused by stiffness, does the strategy convergence noticeably slow down. These results are shown in Table 4 and Table respectively.
Timestep öt
0.01s
0.02s
0.04s
0.08s
Simulation Time Converge
5. 76s
3.28s
2.96s
2.80s
Table 4: The simulation time required for policy convergence on the Push task under different 51.
Contact Coefficient
le
1e
1e
1e
10
Simulation Time to Converge 3.20s
3.28s
3.24s
4.80s
Table The simulation time required for policy convergence on the Push task under different p.
24

---

Published as a conference paper at ICLR 2026
A.6.5 COMPUTATIONAL EFFICIENCY
As described in Section the brute-force eval- uation of Equation 3 exhibits time complex- ity of O(T²) with respect to the number of triangular facets, leading to rapidly increas- ing computational cost as the scene becomes more detailed. In contrast, our method achieves linear complexity in the special case outlined in Appendix A.2, with the only overhead stem- ming from nested optimization and hicrarchical blending, compared to the IPC model (Li al. 2020). To evaluate computational efficiency, compared our method, brute-force compu- tation, and IPC on the Billiards benchmark. We recursively subdivided the mesh of each ball to produce increasingly dense scenes, ranging from 512 to 6208 triangles in total. For each resolution, we performed full trajectory opti- mization to generate Figure With brute-force computation, the average time per frame in- creased significantly from 3.35s to 604.50s. In contrast, our BSH-assisted contact model saw a modest increase from 0.6s to 3.18s. These results empirically confirm the expected quadratic growth of brute-force computation and the near-linear scaling of our method. Interestingly, the IPC model (Li et al. 2020) maintained consistent per-frame cost of approximately 0.5s across all res- olutions. This is attributed to the optimizer requiring fewer Newton iterations mesh resolution increases, offsetting the higher evaluation cost of the contact model. We have also compared our performance using more complex, non-convex geometric shapes. As illustrated Figure [4] replace the cubes in our Gather benchmark with bunny meshes, with 2784 triangles in total. In this benchmark, named Gather-Bunny we compare the average per-frame cost using brute-force computation, our method, and the IPC model, where the cost is 309.62s, 3.89s and 3.14s respec- tively. And in T5) we have designed more challenging articulated manipulation task, named Franka robotic arm push bunny to speci-
Figure 14:
The start (top) and end (bottom) frame of the munny-gather benchmark.
Franka Push, where use cfdegree-of-freedeee fied location. This task involves total of 14343 triangles. We also compare the average per-frame cost using brute-force computation, our method, and the IPC model, with computating cost being > 1000s, 7.83s, and 12s, respectively.
Lo00
0.4
0.2
#Stops
Figure 15: The start (left) and end (middle) frame the Franka-Puse benchmark. The convergency history is shown in the left.
25


Reviews:
[{"review_id": "tbNCu6PxAT", "reviewer": "Reviewer_sWeD", "summary": "This paper presents a novel differentiable contact model for rigid-body simulation that guarantees both physical realism and well-behaved gradients. The authors identify key properties - log-barrier, smoothness, non-prehensility, and non-vanishing gradients as necessary for reliable gradient-based optimization in differentiable physics. They propose a new contact potential formulation that satisfies all these properties and introduce a Bounding Sphere Hierarchy (BSH) mechanism for computational efficiency, reducing the cost of evaluating contact interactions. Extensive experiments across contact-rich tasks demonstrate that this model produces smoother gradients, faster convergence, and improved stability compared to existing baselines.", "questions": "1. What does \"GD\" baseline refer to in Figures 4,6,7, and 8?\n2. Given Figure 5, could the authors quantify how their method scales to mesh complexities in realistic robots (>10k triangles)? Is RL or trajectory optimization still feasible?\n3. The method currently applies to rigid-body systems. Do the authors foresee any principled way to extend the non-vanishing gradient property to deformable or soft-body interactions?\n4. The hierarchical blending (Eq. 8) relies on ϵ > 0 to guarantee theorems. How sensitive are results to this choice in practice?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "1. **Strong theoretical foundation** - the work is grounded in rigorous mathematical analysis, formally defining the necessary conditions for a “well-behaved” contact potential. The authors’ proofs (Appendix A.1–A.3) are precise and significantly strengthen the paper’s technical depth.\n2. **Novel and timely contribution** - differentiable contact modeling remains a bottleneck in differentiable simulation. Introducing a model that satisfies both smoothness and non-vanishing gradient properties is highly relevant, especially as differentiable simulators gain traction in RL and robotics (see Werling et al., 2021; Xu et al., 2022; Ye et al., 2024).\n3. **Thorough experimental validation** - the experiments are well-designed and comprehensive, spanning both manipulation and locomotion tasks. The inclusion of both trivial and random initialization demonstrates robustness.\n4. **High-quality writing and presentation**  -the manuscript is clear, well-organized, and logically develops from theoretical properties to practical implementation.\n5. **Impact potential** - the long-range gradient property directly addresses a key failure mode in current differentiable physics pipelines (gradient vanishing in contact-free configurations), with broad implications for differentiable MPC, co-design, and world-model-based RL.", "weaknesses": "1. **Terminology ambiguity** - the use of “log-barrier” as a property (Property 3.1) is somewhat confusing. In optimization literature, log-barrier refers to a method, not an inherent property. Reframing this as a “barrier potential” or “barrier-form property” would avoid conceptual conflation.\n2. **Referencing and clarity**\n    * Propositions 3.3a and 3.3b are not explicitly labeled in Table 1, which may confuse readers unfamiliar with the notation.\n    * The color scheme in Figure 2 is not self-explanatory; a legend or caption clarification would improve readability.\n    * Figure 3 is visually appealing but does little to elucidate the core algorithmic contribution; it could be condensed or better explained.\n3. **Accessibility for newcomers** - since the proposed formulation heavily builds on Liang et al. (2024) and Ye et al. (2024), a brief, self-contained overview of those works would make this paper more approachable.\n4. **Presentation details** - several figures (especially in Sections 5–6) are too small to read comfortably. As a general guideline, the figure text size should match the main body font.\n5. **Limited domain variety** - while the manipulation and locomotion experiments are strong, an articulated manipulation case (e.g., as in Srinivasan et al., 2024) would better demonstrate real-world robotic applicability.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "q6KP6ADQ3O", "reviewer": "Reviewer_h8Nt", "summary": "Section 2 provides a brief overview on differential simulators in the field of computer graphics. In Section 3, the authors discuss several desirable properties for gradient-based optimization of differentiable simulators. In particular, the authors point out that a contact potential should **approach infinity when the signed distance approaches zero** (aka follow a \"log-barrier\"), be **twice-differentiable** (aka \"smooth\"), **only create pushing forces** (aka \"non-prehensible\"), and **provide non-zero gradients** even if geometries do not touch (\"non-vanishing\"). In **Section 4**, the authors extend the contact potential to fulfill several desirable properties for gradient-based optimization of differentiable simulators. In particular, the authors modify the contact potential function of (Liang et al., 2024; Ye et al., 2024) to have global support. However, computing the proposed potential function\nP is not efficient as it requires computing the contact potential between each pair of\ndisjoint triangles **increasing computation quadratically with the number of triangles**. Therefore, the authors suggest in **Section 5**  to interpolate between two potential function as a function of signed distance. If the objects are close the previously derived potential function is used while if the objects are far apart then a faster to compute potential is evaluated that uses a bounding sphere hierarchy to make  the algorithm computational tractable. In **Section 6** the authors show numerous optimizations of different simulation scenes involving simulator gradients.", "questions": "- Q1: The paper https://github.com/taichi-dev/difftaichi (not cited by the authors) discusses how time-discretization causes gradient errors for stiff contact simulations. How does increasing contact stiffness or integration times affect the gradient quality of the proposed simulator?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "The method developed for smoothly interpolating between different force potentials is quite clever and will definitely inspire follow up work in the field. The derivation of the potential functions appears to be sound.\n\n- Figure 2 is quite nice.\n- Nice experiments (though \"only\" in simulation).", "weaknesses": "- **C1: The notation is a bit convoluted and the writing could be clearer.**\n  Some suggestions:\n\t- You could simply define index $i$ to belong to vertices I and $j$ to vertices J, then you have $f_i$ instead of $f_{i \\in I}^{I \\cup J}$ and $x_i \\in x_{i \\in I}$. Also you could define $\\mathcal{P}^{i \\cup J}$ instead of $\\mathcal{P}$ (as it is clear that the potential function acts between sets of vertices).\n\t- I am not sure what Definition 3.1 is defining and why it required the use of a \"Definition\" statement.\n\t- \"Lemma 4.1\" seems also unnecessary. Also why define properties through numbers (e.g. Property 3.1, 3.2, 3.3, 3.4)?\n     - \"P(x)=\\infty iff x \\in C_obs\" seems imprecise/wrong.\n\t\n- **C2: Adding naive \"global support\" of the contact potential degrades simulation realism.** If you extend the force potential to provide global support then you apply contact forces in the simulation between object that do not touch. While these forces are small, they may notably alter the simulation. For simulation of real-world scene, extending global support could introduce a reality gap and adversely effect control synthesis. That said, in the videos supplied with the supplementary material,  I could not observe such artifacts. That said, the arxiv paper [Hard Contacts with Soft Gradients: Refining Differentiable Simulators for Learning and Control](https://arxiv.org/abs/2506.14186) might be interesting to the authors (this is an arxiv paper so there is no need to discuss it in the author's paper). In principle, the authors could use straight-through estimation to only use the contact potential with global support for the gradient computation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NiMYa0NGk3", "reviewer": "Reviewer_WJTs", "summary": "This paper argues that the choice of contact models for resolving collisions plays a significant role in the quality of gradients computed by differentiable rigid body simulators. They establish necessary properties that a well-behaved contact model should possess, i.e. log-barrier, smoothness, non-prehensile and non-vanishing. The authors then propose a practical contact model which follows all of these properties, and further modify it to improve upon efficiency, while ensuring that the modified contact model is still well-behaved. The proposed method is demonstrated on a variety of gradient based optimization tasks, each of which involves learning a sequence of control signals on controlled objects, which when interacting with other target objects through the proposed contact model, result in a desired motion through simulation. Time complexity analysis is also provided for select tasks.", "questions": "•\tTypo: In Figure 1, the caption should describe a hexagon, while it currently describes a pentagon\n\t•\tIn Section 4, in the expression for computing the contact potential between two triangles, should the RHS be (min L) instead of (argmin L)? \n\t•\tWhat do you use for the value of mu, the coefficient of the contact model P? \n\t•\tWhen the set of vertices are far enough, i..e. at distance d_2 in the notation of the paper, can we simply use the closed form centered potential as in eqn. (5), instead of blending it with the exact potential at d_1? This might be feasible considering the centered potential satisfies all the properties of a well-behaved contact potential as argued in the proof of Corollary 5.2. The general question here is: what is the motivation behind using blending?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "•\tThe theoretical contributions of the paper are quite dense and robust, and they cover all aspects of the proposed methodology in detail, with relevant proofs in the appendix. \n\t•\tThe tasks in the provided demos are quite varied and the corresponding results appear to be convincing.", "weaknesses": "- The exposition can be heavy in notation sometimes, making it hard to understand. \n- Although the authors provide theoretical comparisons with other contact models, it would be even better if the comparison can be extended to the demos in the evaluation section. This might provide valuable insights on the importance of each of the properties described in Section 3.\n- The authors might want to provide an overview of the methodology in terms of an algorithm pseudocode and/or implementation details, for exposition and reproducibility.  \n\nNOTE:   I reviewed an earlier version of this paper.  Most of the text including typos in figures and equations are exactly the same as last time.   I'd like to ask the authors to please at least correct these typos from previous reviews.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SE28zAgKyC", "reviewer": "Reviewer_5Qbn", "summary": "In this paper, the authors formulate four properties a contact model must possess in order to have well-behaved gradients. Those properties are: 1. no interpenetration between contacting bodies, 2. smoothness, 3. unilateral normal forces, 4. non-vanishing gradients. The authors then propose a contact formulation based on potential functions for tedrahedal meshes that satisfies these properties. To increase the computational efficiency of their model, the authors propose to blend between different potentials, i.e., between the exact potential and an approximation based on bounding spheres. The authors evaluate their contact model for different control tasks in simulation and compare to two other contact models and randomized smoothing. The results show that their model provides more informative gradients and enables successful optimization.", "questions": "Thank you for the interesting paper! Maybe I misunderstand something here, but I was confused that the rod overcomes the box to push from the other side in the Push experiment. In my mind, even with contact gradients from a distance, the gradients from the side of the box that is closer to the rod in the beginning should always be stronger than those of the further side. What does the gradient look like if the contact surfaces are \"behind\" each other?", "limitations": "", "rating": 2, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "-The paper rigurously motivates the proposed properties. \n\n-The method (including the performance enhancements) is mathematically described in great detail. \n\n-The results clearly show the gradient quality of the proposed contact model. \n\n-The contact formulation using hierarchical potential blending is novel and computationally efficient.", "weaknesses": "-The claimed contribution of non-vanishing gradients has been studied several times before and is not novel. E.g. see [1,2]\n\n-The baselines the contact model is compared to are limited and not representative of the relevant models in robotics. E.g. see [3,4]\n\n-The supplementary videos should be presented with more effort, one could at least edit them together into a single video.\n\n-My biggest concern with this paper is that physical accuracy is completely neglected. In the beginning, the authors acknowledge that there is a trade off between physical accuracy and well-behaved gradient information. However, they never investigate this trade off and solely focus on gradient quality in the experiments. In my opinion, very simple, smooth contact models can provide the same gradient quality achieved here with much less computation. The advantage of the proposed model could be increased physical realism but without experiments this is hard to judge. The impact and utility of the proposed method are very unclear. \n\n\n[1] Turpin, Dylan, et al. \"Grasp’d: Differentiable contact-rich grasp synthesis for multi-fingered hands.\" European Conference on Computer Vision. Cham: Springer Nature Switzerland, 2022.\n\n[2] Schwarke, Clemens, et al. \"Learning Deployable Locomotion Control via Differentiable Simulation.\" 9th Annual Conference on Robot Learning. 2025.\n\n[3] Todorov, Emanuel, Tom Erez, and Yuval Tassa. \"Mujoco: A physics engine for model-based control.\" 2012 IEEE/RSJ international conference on intelligent robots and systems. IEEE, 2012.\n\n[4] Howell, Taylor A., et al. \"Dojo: A differentiable simulator for robotics.\" arXiv preprint arXiv:2203.00806 9.2 (2022): 4.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)