Write a program whose nontermination is independent of Peano arithmetic Challenge

Write a program P, taking no input, such that the proposition “the execution of P eventually terminates” is independent of Peano arithmetic.

Formal rules

(In case you are a mathematical logician who thinks the above description is too informal.)

One can, in principle, convert some Universal Turing machine U (e.g. your favorite programming language) into a Peano arithmetic formula HALT over the variable p, where HALT(p) encodes the proposition “U terminates on the program (Gödel-encoded by) p”.  The challenge is to find p such that neither HALT(p) nor ¬HALT(p) can be proven in Peano arithmetic.

You may assume your program runs on an ideal machine with unlimited memory and integers/pointers big enough to access it.

Example

To see that such programs exist, one example is a program that exhaustively searches for a Peano arithmetic proof of 0 = 1.  Peano arithmetic proves that this program halts if and only if Peano arithmetic is inconsistent.  Since Peano arithmetic is consistent but cannot prove its own consistency, it cannot decide whether this program halts.

However, there are many other propositions independent of Peano arithmetic on which you could base your program.

Motivation

This challenge was inspired by a new paper by Yedidia and Aaronson (2016) exhibiting a 7,918-state Turing machine whose nontermination is independent of ZFC, a much stronger system than Peano arithmetic.  You may be interested in its citation [22].  For this challenge, of course, you may use your programming language of choice in place of actual Turing machines.
Haskell, 838 bytes

“If you want something done, …”

import Control.Monad.State
data T=V Int|T:$T|A(T-&gt;T)
g=guard
r=runStateT
s!a@(V i)=maybe a id$lookup i s
s!(a:$b)=(s!a):$(s!b)
s@((i,_):_)!A f=A(\a-&gt;((i+1,a):s)!f(V$i+1))
c l=do(m,k)&lt;-(`divMod`sum(1&lt;$l)).pred&lt;$&gt;get;g$m&gt;=0;put m;l!!fromEnum k
i&amp;a=V i:$a
i%t=(:$).(i&amp;)&lt;$&gt;t&lt;*&gt;t
x i=c$[4%x i,5%x i,(6&amp;)&lt;$&gt;x i]++map(pure.V)[7..i-1]
y i=c[A&lt;$&gt;z i,1%y i,(2&amp;)&lt;$&gt;y i,3%x i]
z i=(\a e-&gt;[(i,e)]!a)&lt;$&gt;y(i+1)
(i?h)p=c[g$any(p#i)h,do q&lt;-y i;i?h$q;i?h$1&amp;q:$p,do f&lt;-z i;a&lt;-x i;g$p#i$f a;c[i?h$A f,do b&lt;-x i;i?h$3&amp;b:$a;i?h$f b],case p of A f-&gt;c[(i+1)?h$f$V i,do i?h$f$V 7;(i+1)?(f(V i):h)$f$6&amp;V i];V 1:$q:$r-&gt;c[i?(q:h)$r,i?(2&amp;r:h)$V 2:$q];_-&gt;mzero]
(V a#i)(V b)=a==b
((a:$b)#i)(c:$d)=(a#i)c&amp;&amp;(b#i)d
(A f#i)(A g)=f(V i)#(i+1)$g$V i
(_#_)_=0&lt;0
main=print$(r(8?map fst(r(y 8)=&lt;&lt;[497,8269,56106533,12033,123263749,10049,661072709])$3&amp;V 7:$(6&amp;V 7))=&lt;&lt;[0..])!!0


Explanation

This program directly searches for a Peano arithmetic proof of 0 = 1.  Since PA is consistent, this program never terminates; but since PA cannot prove its own consistency, the nontermination of this program is independent of PA.

T is the type of expressions and propositions:


A P represents the proposition ∀x [P(x)].
(V 1 :$ P) :$ Q represents the proposition P → Q.
V 2 :$ P represents the proposition ¬P.
(V 3 :$ x) :$ y represents the proposition x = y.
(V 4 :$ x) :$ y represents the natural x + y.
(V 5 :$ x) :$ y represents the natural x ⋅ y.
V 6 :$ x represents the natural S(x) = x + 1.
V 7 reprents the natural 0.


In an environment with i free variables, we encode expressions, propositions, and proofs as 2×2 integer matrices [1, 0; a, b], as follows:


M(i, ∀x [P(x)]) = [1, 0; 1, 4] ⋅ M(i, λx [P(x)])
M(i, λx [F(x)]) = M(i + 1, F(x)) where M(j, x) = [1, 0; 5 + i, 4 + j] for all j > i
M(i, P → Q) = [1, 0; 2, 4] ⋅ M(i, P) ⋅ M(i, Q)
M(i, ¬P) = [1, 0; 3, 4] ⋅ M(i, P)
M(i, x = y) = [1, 0; 4, 4] ⋅ M(i, x) ⋅ M(i, y)
M(i, x + y) = [1, 0; 1, 4 + i] ⋅ M(i, x) ⋅ M(i, y)
M(i, x ⋅ y) = [1, 0; 2, 4 + i] ⋅ M(i, x) ⋅ M(i, y)
M(i, S x) = [1, 0; 3, 4 + i] ⋅ M(i, x)
M(i, 0) = [1, 0; 4, 4 + i]
M(i, (Γ, P) ⊢ P) = [1, 0; 1, 4]
M(i, Γ ⊢ P) = [1, 0; 2, 4] ⋅ M(i, Q) ⋅ M(i, Γ ⊢ Q) ⋅ M(i, Γ ⊢ Q → P)
M(i, Γ ⊢ P(x)) = [1, 0; 3, 4] ⋅ M(i, λx [P(x)]) ⋅ M(i, x) ⋅ [1, 0; 1, 2] ⋅ M(i, Γ ⊢ ∀x P(x))
M(i, Γ ⊢ P(x)) = [1, 0; 3, 4] ⋅ M(i, λx [P(x)]) ⋅ M(i, x) ⋅ [1, 0; 2, 2] ⋅ M(i, y) ⋅ M(i, Γ ⊢ y = x) ⋅ M(i, Γ ⊢ P(y))
M(i, Γ ⊢ ∀x, P(x)) = [1, 0; 8, 8] ⋅ M(i, λx [Γ ⊢ P(x)])
M(i, Γ ⊢ ∀x, P(x)) = [1, 0; 12, 8] ⋅ M(i, Γ ⊢ P(0)) ⋅ M(i, λx [(Γ, P(x)) ⊢ P(S(x))])
M(i, Γ ⊢ P → Q) = [1, 0; 8, 8] ⋅ M(i, (Γ, P) ⊢ Q)
M(i, Γ ⊢ P → Q) = [1, 0; 12, 8] ⋅ M(i, (Γ, ¬Q) ⊢ ¬P)


The remaining axioms are encoded numerically and included in the initial environment Γ:


M(0, ∀x [x = x]) = [1, 0; 497, 400]
M(0, ∀x [¬(S(x) = 0)]) = [1, 0; 8269, 8000]
M(0, ∀x ∀y [S(x) = S(y) → x = y]) = [1, 0; 56106533, 47775744]
M(0, ∀x [x + 0 = x]) = [1, 0; 12033, 10000]
M(0, ∀y [x + S(y) = S(x + y)]) = [1, 0; 123263749, 107495424]
M(0, ∀x [x ⋅ 0 = 0]) = [1, 0; 10049, 10000]
M(0, ∀x ∀y [x ⋅ S(y) = x ⋅ y + x]) = [1, 0; 661072709, 644972544]


A proof with matrix [1, 0; a, b] can be checked given only the lower-left corner a (or any other value congruent to a modulo b); the other values are there to enable composition of proofs.

For example, here is a proof that addition is commutative.


M(0, Γ ⊢ ∀x ∀y [x + y = y + x]) = [1, 0; 6651439985424903472274778830412211286042729801174124932726010503641310445578492460637276210966154277204244776748283051731165114392766752978964153601068040044362776324924904132311711526476930755026298356469866717434090029353415862307981531900946916847172554628759434336793920402956876846292776619877110678804972343426850350512203833644, 14010499234317302152403198529613715336094817740448888109376168978138227692104106788277363562889534501599380268163213618740021570705080096139804941973102814335632180523847407060058534443254569282138051511292576687428837652027900127452656255880653718107444964680660904752950049505280000000000000000000000000000000000000000000000000000000]


You can verify it with the program as follows:

*Main&gt; let p = A $ \x -&gt; A $ \y -&gt; V 3 :$ (V 4 :$ x :$ y) :$ (V 4 :$ y :$ x)
*Main&gt; let a = 6651439985424903472274778830412211286042729801174124932726010503641310445578492460637276210966154277204244776748283051731165114392766752978964153601068040044362776324924904132311711526476930755026298356469866717434090029353415862307981531900946916847172554628759434336793920402956876846292776619877110678804972343426850350512203833644
*Main&gt; r(8?map fst(r(y 8)=&lt;&lt;[497,8269,56106533,12033,123263749,10049,661072709])$p)a :: [((),Integer)]
[((),0)]


If the proof were invalid you would get the empty list instead.