For every c>0 the optimized deterministic PBT entanglement-fidelity limit phi(c)=lim_{d->infinity} F_d^*(max(1,floor(c d^2))) exists. It equals c for 0<c<=1/4. For c>1/4, define I_j(mu)=pi^-1 integral_0^pi (cos p-mu)_+^j dp for j=1,2, and let mu be the unique point in (-1,1) satisfying c+1/2=I_2(mu)/(2 I_1(mu)^2). Then phi(c)=[mu+I_2(mu)/I_1(mu)]^2. This addresses the entire literal crossover question; literature novelty and external expert confirmation are unverified.

We use the published finite-resource characterization of optimized deterministic PBT, F_d^*(n)=||D_{d,n}||^2/d^2. Here D_{d,n} maps a partition of n with at most d rows to the sum of its immediate predecessors. Indeed D*D has diagonal entry equal to the number of removable boxes and off-diagonal entry one exactly for a permissible one-box move. Thus this is the teleportation matrix, with the entanglement-fidelity normalization in the original question. All subsequent limits are joint limits with d tending to infinity, not substitutions into fixed-dimensional expansions.

1. Fermions and the bounded commutator. On l^2(N_0) put S|x>=|x-1> for x>=1 and S|0>=0, and X|x>=x|x>. On the d-fold exterior power let T=dGamma(S), V=dGamma(X), and Q=V-d(d-1)/2. Under x_1<...<x_d <-> lambda_i=x_{d+1-i}-(d-i), Q is the Young-diagram size. A permissible one-step particle lowering has positive coefficient: it never crosses a particle. Thus T on grade Q=n is exactly D_{d,n}, and T annihilates grade zero. Also ||T||<=d and
[T,T*]=dGamma(SS*-S*S)=dGamma(|0><0|)=:K,
where 0<=K<=I because one fermionic site has occupation zero or one. These are identities on the full wedge; Q and V are unbounded, but every vector used initially below has finite support.

2. Transfer from a fixed grade to an expectation. Write A=T*T and let psi be a unit top eigenvector of A in grade n, with eigenvalue a^2=||D_{d,n}||^2. The grade is finite dimensional and invariant under A. Fix m independently of d and put v_j=(T*)^j psi for 0<=j<=m. Since [A,T*]=T*K, telescoping gives
||(A-a^2)v_j|| <= j d^j.
Suppose along a subsequence a/d is bounded below by some eta>0. Induction, using
||v_{j+1}||^2=a^2||v_j||^2+<v_j,(A-a^2)v_j>+<v_j,Kv_j>,
shows ||v_j||=a^j(1+O_{m,eta}(d^-2)). To see the error explicitly, if ||v_j|| is of order d^j then the last two terms have absolute value at most j d^j||v_j||+||v_j||^2=O_{m,eta}(d^{2j}), while the first is of order d^{2j+2}. This starts with ||v_0||=1 and establishes the induction, including nonvanishing.
The normalized w_j=v_j/||v_j|| lie in orthogonal grades n+j, and
<w_j,Tw_{j+1}>=||v_{j+1}||/||v_j||=a(1+O_{m,eta}(d^-2)).
Set alpha_j=sqrt(2/(m+2)) sin(pi(j+1)/(m+2)), and z=sum_{j=0}^m alpha_j w_j. The elementary path-matrix sine identity yields
<z,Re(T)z>/d=(a/d)cos(pi/(m+2))+O_{m,eta}(d^-2),
and <z,Qz><=n+m. If instead a/d tends to zero, every nonnegative limiting upper bound below already holds. Consequently an asymptotic upper bound on Re(T)/d for states of mean grade at most n+m bounds the fixed-grade singular value after first d->infinity and then m->infinity. This is a proved transfer, not an assumption of ensemble equivalence.

3. An elementary one-particle trace limit. For fixed b>0 and real mu let h_{b,d}=Re(S)-bX/d. It is self-adjoint on the domain of X, is bounded above, and has compact resolvent. We claim
lim_{d->infinity} d^-1 Tr(h_{b,d}-mu)_+ = I_2(mu)/(2b),
where I_j(mu)=pi^-1 integral_0^pi (cos p-mu)_+^j dp.
Here is a direct proof. Choose C>max(0,(1-mu)/b)+1 and remove the hopping edge between the sites below floor(Cd) and the remaining half-line. The difference has trace norm one. On the tail, Re(S)<=I and X/d>=floor(Cd)/d, so the tail operator is strictly below mu for all sufficiently large d and contributes no positive trace. For self-adjoint operators bounded above with finite positive trace, the variational characterization Tr(B)_+=sup_{0<=P<=I, P finite rank}Tr(PB) shows that a trace-class perturbation changes this trace by at most its trace norm. Thus the removed edge costs at most one.
Partition the retained interval into a fixed number M of blocks with endpoints approximating a partition of [0,C]. Removing the M-1 internal hopping edges costs at most M-1. On each block, replace X/d by either endpoint; the two resulting operators order the original block. The gap in potential is at most b times the mesh plus O(1/d). The corresponding normalized positive traces differ by at most C b times the mesh plus o(1), since positive-part traces are monotone under operator order and a scalar shift epsilon changes the trace in dimension L by at most L|epsilon|. A length-L block with constant potential -by has eigenvalues cos(k pi/(L+1))-by, k=1,...,L. First let d grow with this partition fixed; these sums are ordinary Riemann sums. Then let the mesh tend to zero. This gives
(1/pi) integral_0^C integral_0^pi (cos p-by-mu)_+ dp dy.
The integrand vanishes beyond C, and integrating y yields I_2(mu)/(2b), as claimed. Every cutting error was O(M)/d before refinement; no uniform fixed-d asymptotic is being invoked.

4. Universal upper bound. For any normalized d-fermion state z of finite mean V, the fermionic variational principle bounds the expectation of dGamma(h_{b,d}) by the sum of the d largest eigenvalues of h_{b,d}. This also follows from the one-particle density matrix 0<=gamma<=I, Tr gamma=d, by filling the d largest levels. For any mu, their sum is at most d mu+Tr(h_{b,d}-mu)_+. Therefore
<Re(T)>/d <= mu+d^-1 Tr(h_{b,d}-mu)_+ + b <V>/d^2.
For n=N_d=max(1,floor(cd^2)), the states z from step 2 have <V>/d^2<=c+1/2+o(1). Along every positive subsequential limit of a/d, steps 2 and 3, followed by m->infinity, give
limsup a/d <= mu+b(c+1/2)+I_2(mu)/(2b)
for each fixed b>0 and mu. The same assertion is trivial on subsequences with limit zero whenever the chosen right side is nonnegative.

5. The parameters and the resulting upper bound. Write M=c+1/2. If 0<c<=1/4 choose mu=-1/(2sqrt(c))<=-1 and b=I_1(mu)=-mu. Since I_2(mu)=mu^2+1/2, M=I_2/(2I_1^2), and the preceding bound equals
q(c)=mu+I_2(mu)/I_1(mu)=sqrt(c).
If c>1/4 choose the unique mu in (-1,1) with M=I_2(mu)/(2I_1(mu)^2), and again b=I_1(mu). This gives the bound q(c)=mu+I_2/I_1. Existence, uniqueness and continuity require no variational guess: for -1<mu<1 put theta=arccos(mu), J=theta/pi. Differentiation gives I_1'=-J and I_2'=-2I_1. Hence the derivative of I_2/(2I_1^2) is (I_2 J-I_1^2)/I_1^3>0 by strict Cauchy-Schwarz on [0,theta]. Its limit at mu=-1 is 3/4. As theta decreases to zero, I_1~theta^3/(3pi), I_2~2theta^5/(15pi), and the ratio is asymptotic to 3pi/(5theta), hence diverges. The displayed asymptotics follow by Taylor expansion of cosine on this shrinking interval with uniform remainder. At c=1/4 both definitions give q=1/2. Also q>0 because it equals the positive integral of cos p(cos p-mu)_+ divided by pi I_1; when mu<0, pairing p with pi-p verifies positivity, and when mu>=0 positivity is immediate. Thus limsup F_d^*(N_d)<=q(c)^2 for every c>0.

6. Matching, explicitly discretizable Slater states. Fix 0<c'<c and choose its parameters mu,b as in step 5. On y>=0 define
f(y)=(1/pi) arccos(mu+by) if -1<mu+by<1,
f(y)=1 if mu+by<=-1, and f(y)=0 if mu+by>=1.
It is continuous and supported on [0,Y], Y=(1-mu)/b. Fubini gives
integral f(y)dy=I_1/b=1,
integral y f(y)dy=I_2/(2b^2)=c'+1/2,
integral sin(pi f(y))/pi dy=Icos/b=q(c'),
where Icos=pi^-1 integral cos p(cos p-mu)_+ dp=I_2+mu I_1. The last identity can also be seen by integrating the kinetic energy cos p over the occupied region 0<=p<=pi f(y).

Choose a finite partition 0=y_0<...<y_M=Y and replace f by its cell averages f_j. Then sum (y_j-y_{j-1})f_j=1 exactly. Its cell-midpoint first moment tends to integral y f, and sum (y_j-y_{j-1})sin(pi f_j)/pi tends to q(c') as the mesh vanishes, by continuity. Choose the mesh small enough that its first moment minus 1/2 is strictly less than c.
For large integer d, form consecutive disjoint site blocks with endpoints floor(d y_j) and lengths L_j. Put k_j fermions in block j, using its first k_j normalized sine orbitals. Choose integers 0<=k_j<=L_j with sum k_j=d and k_j/d->(y_j-y_{j-1})f_j. Such integers are obtained by rounding and O(M) corrections: total desired count is d+O(M), and at least one positive-width cell has 0<f_j<1, supplying the O(M) spare occupied and vacant positions for all sufficiently large d. Take the Slater determinant of all d orbitals. Different blocks have zero cross-block one-body density, so hopping across their boundaries has zero expectation. The sine formulas give
<T>/d -> sum_j (y_j-y_{j-1})sin(pi f_j)/pi,
<Q>/d^2 -> sum_j (y_j-y_{j-1})(y_{j-1}+y_j)f_j/2 -1/2 <c.
All occupied sites are below Yd+O(1). If P is the one-particle orbital projection, the Slater occupation identities give
Var(Q)=Tr(PX(I-P)X)<=Tr(PX^2)<=d(Yd+O(1))^2=O(d^3).
Thus Pr(Q>N_d)=O(1/d) by the strictly positive mean-size margin and Chebyshev.

7. Transfer the lower bound to exactly N_d ports. Decompose this Slater state psi=sum psi_n by Q. Different T psi_n occupy distinct grades. The operational fidelity F_d^*(n) is nondecreasing in n: an extra port may be ignored using zero measurement effects and a product extension of the resource. Also ||T||<=d. Therefore
|<psi,Tpsi>|^2/d^2 <= ||Tpsi||^2/d^2
 <= F_d^*(N_d) sum_{n<=N_d}||psi_n||^2 + sum_{n>N_d}||psi_n||^2
 <= F_d^*(N_d)+Pr(Q>N_d).
The grade-zero term vanishes. First d->infinity gives the squared block-kinetic lower bound. Then refine the finite partition; it remains inside the positive margin for all sufficiently fine partitions, giving liminf F_d^*(N_d)>=q(c')^2. Finally c' increases to c; continuity from step 5 yields liminf>=q(c)^2. Combined with the upper bound this proves existence and the claimed value for every c>0. All mesh, ladder and margin limits are taken only after the growing-d limit.

For convenient evaluation, if -1<mu<1 and theta=arccos(mu), then
I_1=(sin(theta)-mu theta)/pi,
I_2=((1/2+mu^2)theta-(3/2)mu sin(theta))/pi.
These elementary expressions specify the answer without a numerical fitting step. The supplementary finite checker only tests the exact branching commutator in small dimensions; it does not machine-certify the asymptotic argument, whose proof is above.
