def realize(E):
"""one-mode realization phi_i = cu_i q + cw_i p of a rank-2 antisymmetric E"""
T, Z = schur(E, output='real')
k = next(k for k in range(len(E) - 1) if abs(T[k, k + 1]) > 1e-12)
b, uu, ww = T[k, k + 1], Z[:, k], Z[:, k + 1]
if b < 0:
b, uu, ww = -b, ww, uu
cu, cw = np.sqrt(b) * uu, np.sqrt(b) * ww
assert np.abs(np.outer(cu, cw) - np.outer(cw, cu) - E).max() < 1e-12
return cu, cw
ncut = 140
a = ladder(ncut)
qo, po = (a + a.T) / np.sqrt(2), 1j * (a.T - a) / np.sqrt(2)
edges = [-np.inf, -0.7, 0.0, 0.7, np.inf]
def run_order(order, eigs, mult):
"""apply the maps rho -> V (M o V* rho V) V* in the given order to the vacuum; trace"""
rho = np.zeros((ncut, ncut), complex)
rho[0, 0] = 1
for m in order:
w_, V_ = eigs[m]
rho = V_ @ ((V_.conj().T @ rho @ V_) * mult(m, w_)) @ V_.conj().T
return np.trace(rho)
def bin_mult(bins, sg):
def mb(m, w_):
dd, mm = w_[:, None] - w_[None, :], (w_[:, None] + w_[None, :]) / 2
lo, hi = edges[bins[m]], edges[bins[m] + 1]
mass_in_bin = Phi((hi - mm) / sg[m]) - Phi((lo - mm) / sg[m])
return np.exp(-dd**2 / (8 * sg[m]**2)) * mass_in_bin
return mb
def bin_probability(rho, eig, b, sgy):
"""probability that a Gaussian collapse of width sgy lands in bin b (eig: its observable)"""
w_, V_ = eig
dd, mm = w_[:, None] - w_[None, :], (w_[:, None] + w_[None, :]) / 2
in_bin = Phi((edges[b + 1] - mm) / sgy) - Phi((edges[b] - mm) / sgy)
M = np.exp(-dd**2 / (8 * sgy**2)) * in_bin
return np.trace((V_.conj().T @ rho @ V_) * M).real
def binned_law(order, eigs, sg, k):
law = np.zeros((4,) * k)
for bins in np.ndindex(*law.shape):
law[bins] = run_order(order, eigs, bin_mult(bins, sg)).real
return law
found, cf_dev, ord_dev, norm_dev = 0, 0.0, 0.0, 0.0
while found < 4:
pts = rng.uniform(0, 1.6, size=(3, 2))
if topo(relation(pts, dl)) is None:
continue
R, E = geometry(pts, dl, scale)
orders = list({tuple(topo(R, rng)) for _ in range(30)})
if np.sum(np.abs(E) > 1e-3) < 4 or len(orders) < 2:
continue
found += 1
cu, cw = realize(E)
eigs = [eigh(cu[i] * qo + cw[i] * po) for i in range(3)]
Sig = 0.5 * np.outer(cu, cu) + 0.5 * np.outer(cw, cw)
sg = np.array([0.8, 0.9, 1.0])
for s in rng.normal(size=(4, 3)):
def ms(m, w_, s=s):
dd, mm = w_[:, None] - w_[None, :], (w_[:, None] + w_[None, :]) / 2
return np.exp(-dd**2 / (8 * sg[m]**2) + 1j * s[m] * mm - s[m]**2 * sg[m]**2 / 2)
for o in orders:
closed = np.exp(-0.5 * s @ cov_closed(list(o), E, Sig, sg) @ s)
cf_dev = max(cf_dev, abs(run_order(o, eigs, ms) - closed))
laws = [binned_law(o, eigs, sg, 3) for o in orders]
norm_dev = max(norm_dev, max(abs(l.sum() - 1) for l in laws))
ord_dev = max(ord_dev, max(np.abs(l - laws[0]).max() for l in laws))
print(' (d) 4 three-collapse configurations, one-mode realization, Fock cutoff %d:' % ncut)
print(' characteristic function, direct vs closed form: max dev %.1e' % cf_dev)
print(' 64-cell binned history laws: max difference between linear extensions %.1e,'
% ord_dev)
print(' normalization error %.1e' % norm_dev)
pts = np.array([[0.3, 0.1], [0.55, 0.25]]) # two overlapping diamonds
R, E = geometry(pts, dl, scale)
cu, cw = realize(E)
eigs = [eigh(cu[i] * qo + cw[i] * po) for i in range(2)]
sg = np.array([0.8, 0.9])
rho = np.zeros((ncut, ncut), complex)
rho[0, 0] = 1
A2 = cu[1] * qo + cw[1] * po
rho1 = dephase(rho, cu[0] * qo + cw[0] * po, sg[0])
var_after = (np.trace(rho1 @ A2 @ A2) - np.trace(rho1 @ A2)**2).real
var_alone = (np.trace(rho @ A2 @ A2) - np.trace(rho @ A2)**2).real
print(' (e) overlapping pair, E(h1, h2) = %.6f:' % E[0, 1])
print(' increase of Var(y2) caused by collapse 1: %.10f; E^2/(4 sigma1^2) = %.10f'
% (var_after - var_alone, E[0, 1]**2 / (4 * sg[0]**2)))
l12, l21 = binned_law((0, 1), eigs, sg, 2), binned_law((1, 0), eigs, sg, 2)
print(' binned joint laws for the two orders: total-variation distance %.4f'
% (0.5 * np.abs(l12 - l21).sum()))
A1 = cu[0] * qo + cw[0] * po
r12 = dephase(dephase(rho, A1, sg[0]), A2, sg[1])
r21 = dephase(dephase(rho, A2, sg[1]), A1, sg[0])
print(' unread channels in the two orders: max difference of the final states %.1e'
% np.abs(r12 - r21).max())
eta, nbad, nchg, maxd = 0.6, 0, 0, 0.0
while nbad < 300:
pts = rng.uniform(0, 2.0, size=(rng.poisson(6.0), 2))
if len(pts) < 2 or topo(relation(pts, dl)) is not None:
continue
n = len(pts)
R, E = geometry(pts, dl, scale)
sg = np.full(n, 0.8)
Sig = 0.5 * absE(E) + 0.3 * np.eye(n)
o1 = list(np.argsort(pts[:, 0]))
o2 = list(np.argsort(np.cosh(eta) * pts[:, 0] - np.sinh(eta) * pts[:, 1]))
dC = np.abs(cov_closed(o1, E, Sig, sg) - cov_closed(o2, E, Sig, sg)).max()
nbad, nchg, maxd = nbad + 1, nchg + (dC > 1e-12), max(maxd, dC)
print(' (f) %d non-orderable configurations ordered by centre time in two frames' % nbad)
print(' (relative rapidity 0.6): laws differ in %d, max covariance difference %.3f'
% (nchg, maxd))
print(' (g) centre-time order is a linear extension of the causal order in %d of the %d'
% (t_ext, nconf))
print(' orderable configurations of item (c)')
# (h) cosmological-time-ordered collapses on a non-orderable Sorkin-type configuration:
# j (tested operation), c1 and c2 (overlapping law collapses), y (target), applied in order
# of t. Light-cone centres (u, v) give j < c1 < y causally, j and y causally separated.
uvc = np.array([[0.0, 0.0], [-0.525, 1.2], [-0.225, 1.05], [-1.05, 2.3]])
pts = np.column_stack([(uvc[:, 0] + uvc[:, 1]) / 2, (uvc[:, 1] - uvc[:, 0]) / 2])
R, E = geometry(pts, dl, 8.0)
print(' (h) order j, c1, c2, y by t; c1 and c2 overlap; E(h_j, h_c1) = %.4f, E(h_c1, h_y) = %.4f,'
% (E[0, 1], E[1, 3]))
print(' E(h_c1, h_c2) = %.4f, E(h_j, h_y) = %.1f' % (E[1, 2], abs(E[0, 3])))
gl = lambda x: (1 + np.tanh(2 * x)) / 2
for nc in (36, 44):
a1, a2 = modes(nc, 2)
q1, q2 = [(a + a.T) / np.sqrt(2) for a in (a1, a2)]
p1, p2 = [1j * (a.T - a) / np.sqrt(2) for a in (a1, a2)]
# two-mode realization of the commutation data with phi_j on mode 1 and phi_y on mode 2
phis = [q1, E[1, 2] / E[0, 2] * q1 + E[0, 1] * p1 - E[1, 3] * p2,
E[0, 2] * p1 - E[2, 3] * p2, q2]
for i in range(4):
for k in range(4):
c = (phis[i] @ phis[k] - phis[k] @ phis[i])[:8, :8]
assert np.abs(c - 1j * E[i, k] * np.eye(8)).max() < 1e-9
ey = eigh(phis[3])
w1, V1 = eigh(phis[1])
w0, V0 = eigh(phis[0])
def target_law(at_j, luders):
rho = np.zeros((nc * nc, nc * nc), complex)
rho[0, 0] = 1
if at_j == 'kick':
U = expm(3j * phis[0])
rho = U @ rho @ U.conj().T
elif at_j == 'collapse':
rho = dephase(rho, phis[0], 0.8)
elif at_j == 'luders': # an operation outside the Weyl-diagonal class
K = [V0 @ np.diag(np.sqrt(gl(w0))) @ V0.conj().T,
V0 @ np.diag(np.sqrt(1 - gl(w0))) @ V0.conj().T]
rho = sum(k @ rho @ k.conj().T for k in K)
if luders:
K = [V1 @ np.diag(np.sqrt(gl(w1))) @ V1.conj().T,
V1 @ np.diag(np.sqrt(1 - gl(w1))) @ V1.conj().T]
rho = sum(k @ rho @ k.conj().T for k in K)
else:
rho = dephase(rho, phis[1], 0.8)
rho = dephase(rho, phis[2], 0.9)
return np.array([bin_probability(rho, ey, b, 1.0) for b in range(4)])
base = target_law(None, False)
dk = np.abs(target_law('kick', False) - base).max()
dc = np.abs(target_law('collapse', False) - base).max()
lb = target_law(None, True)
dlk = np.abs(target_law('kick', True) - lb).max()
print(' cutoff %d: Gaussian c1: change of the target law by a kick at j %.1e,' % (nc, dk))
dlj = np.abs(target_law('luders', False) - base).max()
print(' by a collapse at j %.1e; target law sums to %.12f' % (dc, base.sum()))
print(' unsharp Lueders instrument at j (not Weyl-diagonal): change %.1e' % dlj)
print(' unsharp Lueders c1 (negative control): change by the kick, largest bin %.7f'
% dlk)
print('C4. Gapless pointer realization of the collapses (Proposition 13.31)')
# Each collapse is realized by a pointer (q_k, p_k) with no Hamiltonian, coupled through the
# density p_k h_k(x) phi(x). The matrix M below is Gamma_kl = E_ret(h_k, h_l), the retarded
# part of E, with E_ret(t, x) = -theta(u) theta(v)/2, so that E = Gamma - Gamma^T.
def Mret(pts, dl, scale):
u, v = pts[:, 0] - pts[:, 1], pts[:, 0] + pts[:, 1]
n = len(pts)
M = np.zeros((n, n))
for k in range(n):
for l in range(n):
M[k, l] = -scale**2 / 32 * ((1 + Sfun((u[k] - u[l]) / dl))
* (1 + Sfun((v[k] - v[l]) / dl)))
return M
def ordered_sums(pts, dl, scale, nn, eta):
"""sum of E(x_a - x_b) w_a w_b over node pairs with x_a later than x_b in the time
cosh(eta) t - sinh(eta) x: the phase of the time-ordered product of the couplings"""
z, wz = np.polynomial.legendre.leggauss(nn)
Zu, Zv = np.meshgrid(z, z, indexing='ij')
w = (scale / 2) * np.outer(beta(z) * wz, beta(z) * wz).ravel()
nodes = [((t - x) + dl * Zu.ravel(), (t + x) + dl * Zv.ravel()) for t, x in pts]
T = np.zeros((len(pts), len(pts)))
for k in range(len(pts)):
for l in range(len(pts)):
for c in range(0, len(w), 512):
du = nodes[k][0][c:c + 512, None] - nodes[l][0][None, :]
dv = nodes[k][1][c:c + 512, None] - nodes[l][1][None, :]
later = np.cosh(eta) * (du + dv) - np.sinh(eta) * (dv - du) > 0
Ev = np.where(later, -(np.sign(du) + np.sign(dv)) / 4, 0.0)
T[k, l] += w[c:c + 512] @ Ev @ w
return T