A problem in indexing hub selection for a many-body state problem and erdos number
I'm working on reproducing figure 8 of “Dynamical decoherence of a qubit coupled to a quantum dot or the SYK black hole,”.
The figure is:
In the paper "hub=0 corresponding to the many-body state where first 7 out of 16 orbitals are occupied." The Hamiltonian describing these states is:
I tried a python code to reproduce this figure. This code is very long so I will reproduce here only the part that is convenient for this question, but I can with no problem provide everything if needed. The code is:
"""
import numpy as np
import matplotlib.pyplot as plt
import time
import itertools
from scipy.special import comb
from collections import deque
from itertools import combinations
from matplotlib.ticker import LogLocator
class TBRIMHamiltonian:
"""
Implementation of the spinless TBRIM Hamiltonian following the paper
H = sum_k eps_k c_k^† c_k
+ (4 / sqrt(2 M^3)) sum_{i> k) & 1]
def _generate_Vijkl(self):
"""
Generate antisymmetrized two-body tensor
INCLUDING diagonal interaction terms V_{ij,ij}
"""
V = {}
for i in range(self.M):
for j in range(i + 1, self.M):
# diagonal interaction (density-density)
x_diag = np.random.normal(0.0, self.J)
V[(i, j, i, j)] = x_diag
V[(j, i, j, i)] = x_diag
# off-diagonal scattering
for k in range(self.M):
for l in range(k + 1, self.M):
if (i, j) == (k, l):
continue
x = np.random.normal(0.0, self.J)
V[(i, j, k, l)] = x
V[(j, i, k, l)] = -x
V[(i, j, l, k)] = -x
V[(j, i, l, k)] = x
return V
# Fermionic sign
def _fermion_sign(self, state, annihilate, create):
"""
Compute fermionic sign for applying:
c_{create[0]}^† c_{create[1]}^† c_{annihilate[1]} c_{annihilate[0]}
"""
sign = 1
s = state
# annihilation operators (right to left)
for a in annihilate:
sign *= (-1) ** (bin(s & ((1 << a) - 1)).count("1"))
s &= ~(1 << a)
# creation operators (left to right)
for c in create:
sign *= (-1) ** (bin(s & ((1 << c) - 1)).count("1"))
s |= (1 << c)
return sign
# Matrix element
def matrix_element(self, state_i, state_j):
"""Compute ⟨i|H|j⟩"""
# Diagonal
if state_i == state_j:
occ = self._occupied_orbitals(state_i)
# one-body term
E = sum(self.eps_k[k] for k in occ)
# diagonal interaction term
for i, j in combinations(occ, 2):
if (i, j, i, j) in self.Vijkl:
E += self.prefactor * self.Vijkl[(i, j, i, j)]
return E
# Off-diagonal: must differ by exactly 4 orbitals
diff = state_i ^ state_j
diff_orbs = [k for k in range(self.M) if (diff >> k) & 1]
if len(diff_orbs) != 4:
return 0.0
occ_i = set(self._occupied_orbitals(state_i))
occ_j = set(self._occupied_orbitals(state_j))
annihilate = sorted([k for k in diff_orbs if k in occ_i])
create = sorted([k for k in diff_orbs if k in occ_j])
if len(annihilate) != 2 or len(create) != 2:
return 0.0
i, j = create
k, l = annihilate
if (i, j, k, l) not in self.Vijkl:
return 0.0
sign = self._fermion_sign(state_i, annihilate, create)
return self.prefactor * sign * self.Vijkl[(i, j, k, l)]
# Build full Hamiltonian
def build_full_hamiltonian(self):
"""Construct the full many-body Hamiltonian matrix."""
print("Building Hamiltonian matrix...")
H = np.zeros((self.d, self.d))
for a, state_a in enumerate(self.basis):
H[a, a] = self.matrix_element(state_a, state_a)
for b in range(a + 1, self.d):
val = self.matrix_element(state_a, self.basis[b])
if abs(val) > 1e-14:
H[a, b] = val
H[b, a] = val
return H
class Network:
def build_network(self, H, C):
print(f"Building network with parameter C = {C}...")
d = H.shape[0]
adj = np.zeros((d, d), dtype=bool)
diag = np.diag(H)
# many-body level spacing estimate
bandwidth = np.max(diag) - np.min(diag)
Delta = bandwidth / d if bandwidth > 0 else 1.0 / d
connections = 0
for i in range(d):
for j in range(i+1, d):
denom = max(abs(diag[i] - diag[j]), Delta)
if abs(H[i,j]) > C * denom:
adj[i,j] = True
adj[j,i] = True
connections += 1
print(f"Created {connections} connections "
f"(density : {100*connections/(d*(d-1)/2):.4f}%)")
return adj, diag
def compute_link_distribution(self, adj):
"""Compute the distribution of link numbers"""
link_counts = np.sum(adj, axis=1)
# Statistics
mean_links = np.mean(link_counts)
std_links = np.std(link_counts)
# Histogram
max_links = int(np.max(link_counts))
bins = np.arange(0, max_links+2)
hist, _ = np.histogram(link_counts, bins=bins)
return link_counts, hist, bins, mean_links, std_links
def compute_erdos_numbers(self, adj, hub_index):
"""
Compute Erdös numbers using BFS
(breadth-first search)
"""
n = adj.shape[0]
erdos = -np.ones(n, dtype=int) # -1 means unreachable
erdos[hub_index] = 0
queue = deque([hub_index])
while queue:
current = queue.popleft()
current_dist = erdos[current]
# Find neighbors
neighbors = np.where(adj[current])[0]
for neighbor in neighbors:
if erdos[neighbor] == -1:
erdos[neighbor] = current_dist+1
queue.append(neighbor)
# Statistics
reachable_mask = erdos != -1
reachable_erdos = erdos[reachable_mask]
if len(reachable_erdos) > 1:
mean_erdos = np.mean(reachable_erdos[reachable_erdos > 0])
std_erdos = np.std(reachable_erdos[reachable_erdos > 0])
else:
mean_erdos = 0
std_erdos = 0
reachable_fraction = np.sum(reachable_mask) / n
return erdos, reachable_fraction, mean_erdos, std_erdos
class PlotFigure:
def plot_figure(self, results, param_sets, C_values):
"""
Create all the plots
"""
fig, axes = plt.subplots(2, 2, figsize=(14,10))
# Plot for each C value
for c_idx, C in enumerate(C_values):
ax_erdos = axes[0, c_idx]
ax_links = axes[1, c_idx]
ax_erdos.set_title(f"Erdös Distribution (C={C})", fontsize=12)
ax_erdos.set_xlabel("Erdős number $N_E$")
ax_erdos.set_ylabel("Frequency $N_f(N_E)$")
ax_links.set_title(f"Link Distribution (C={C})", fontsize=12)
ax_links.set_xlabel("Link number $N_l$")
ax_links.set_ylabel("Frequency $N_f(N_l)$")
ax_erdos.set_yscale("log")
ax_links.set_yscale("log")
ax_links.set_xscale("log")
#ax_erdos.set_ylim(bottom=1e-2)
ax_erdos.yaxis.set_major_locator(LogLocator(base=10.0))
ax_erdos.yaxis.set_minor_locator(LogLocator(base=10.0, subs="auto"))
ax_links.yaxis.set_major_locator(LogLocator(base=10.0))
ax_links.yaxis.set_minor_locator(LogLocator(base=10.0, subs="auto"))
ax_links.xaxis.set_major_locator(LogLocator(base=10.0))
ax_links.xaxis.set_minor_locator(LogLocator(base=10.0, subs="auto"))
for params in param_sets:
label = params['label']
color = params['color']
data = results[label][C]
# Plot Erdös distribution
erdos_counts = np.array(data["all_erdos"])
unique_erdos, counts = np.unique(erdos_counts, return_counts=True)
# normalize to probability
counts = counts / n_realizations
# Remove zero counts (log-scale safety)
mask = counts > 0
unique_erdos = unique_erdos[mask]
counts = counts[mask]
if len(unique_erdos) > 0:
ax_erdos.plot(
unique_erdos,
counts,
'o-',
color=color,
#label=f"{label}\nμ={data['mean_erdos']:.2f}±{data['std_erdos']:.2f}"
label=f"{label}\nV={params["V"]:.4f}, J={params["J"]:.4f}"
)
# Plot links distribution
link_counts = np.array(data["all_links"])
# Normalize by maximum possible links
#link_density = link_counts / (len(link_counts) - 1)
# Optional: add tiny jitter to avoid single-point collapse
#jitter = 0.01
jitter = 0.5
#link_density_jittered = link_density + np.random.uniform(-jitter, jitter, size=link_density.shape)
link_counts_jittered = link_counts + np.random.uniform(-jitter, jitter, size=link_counts.shape)
#unique_links, link_hist = np.unique(np.round(link_counts_jittered, 2), return_counts=True)
max_links = int(np.max(link_counts))
bins = np.arange(0, max_links+2)
hist, _ = np.histogram(np.round(link_counts_jittered), bins=bins)
unique_links = bins[:-1]
#link_hist = hist
link_hist = hist.astype(float)
link_hist = link_hist / n_realizations
if len(unique_links) > 0:
ax_links.plot(
unique_links,
link_hist,
's--',
color=color,
#label=f'{label}\nμ={data["mean_links"]:.2f}±{data["std_links"]:.2f}'
label=f"{label}\nV={params["V"]:.4f}, J={params["J"]:.4f}"
)
"""
# Plot links distribution
link_counts =data['link_counts']
unique_links, link_hist = np.unique(link_counts, return_counts=True)
if len(unique_links) > 0:
ax_links.plot(unique_links, link_hist, 's--', color=color,
label=f'{label}\nμ={data["mean_links"]:.2f}±{data["std_links"]:.2f}')
"""
ax_erdos.legend(fontsize=8, loc='best')
ax_links.legend(fontsize=8, loc='best')
ax_erdos.grid(True, alpha=0.3)
ax_links.grid(True, alpha=0.3)
ymin, ymax = ax_erdos.get_ylim()
ax_erdos.set_ylim(1e-2, ymax)
ymin, ymax = ax_links.get_ylim()
ax_links.set_ylim(1e-2, ymax)
plt.suptitle(f"TBRIM Network Analysis, M={M} and L={L} for {n_realizations} realizations", fontsize=14)
plt.tight_layout()
plt.show()
# Small-world properties
print("\n" + "="*60)
print("SMALL-WORLD ANALYSIS")
print("="*60)
for params in param_sets:
label = params['label']
print(f"\n{label}:")
for C in C_values:
data = results[label][C]
adj = data['adj']
n = adj.shape[0]
# Average Path length
sampled_paths = []
for i in range(min(20, n)):
distances = self.bfs_from_node(adj, i)
valid_dist = distances[distances != -1]
if len(valid_dist) > 1:
sampled_paths.append(np.mean(valid_dist[valid_dist > 0]))
avg_L = np.mean(sampled_paths) if sampled_paths else float('inf')
# Connection probability
p = np.mean(adj)
# Random graph comparison
if p > 0:
L_random = np.log(n) / np.log(p*n)
small_world_ratio = avg_L / L_random if L_random > 0 else float('inf')
print(f" C={C}: L={avg_L:.2f}, L_random={L_random:.2f}, ratio={small_world_ratio:.2f}")
if 0.8 < small_world_ratio < 1.2 and data['mean_erdos'] < 4:
print(f" ✓ Small-world property detected!")
else:
print(f" C={C}: No connections")
def bfs_from_node(self, adj, start):
"""BFS to compute distances from start node"""
n = adj.shape[0]
distances = -np.ones(n, dtype=int)
distances[start] = 0
queue = deque([start])
while queue:
current = queue.popleft()
current_dist = distances[current]
neighbors = np.where(adj[current])[0]
for neighbor in neighbors:
if distances[neighbor] == -1:
distances[neighbor] = current_dist+1
queue.append(neighbor)
return distances
class RunModelAndAnalysis:
def run_model(self, M, L, C_values, n_realizations):
"""Run the model :
Build the Hamiltonian
"""
d = int(comb(M, L, exact=True))
print(f"System: M={M}, L={L}, d={d}")
# Define parameter sets
param_sets = [
{"V": 0.0, "J": 1.0, "label": "SYK", "color": "red"},
{"V": np.sqrt(14.0), "J": 0.025, "label": "Weak", "color": "green"},
{"V": np.sqrt(14.0), "J": 0.25, "label": "Moderate", "color": "blue"},
{"V": np.sqrt(14.0), "J": 1.0, "label": "Strong", "color": "purple"}
]
results = {}
for params in param_sets:
print(f"\n{'='*40}")
print(f"Analyzing {params['label']}: V={params['V']}, J={params['J']}")
print(f"\n{'='*40}")
label = params["label"]
results[label] = {
C: {
"all_links": [],
"all_erdos": [],
"reach_frac": []
} for C in C_values
}
for r in range(n_realizations):
print(f" Realization {r+1}/{n_realizations}")
tbrim = TBRIMHamiltonian(
M, L, params['V'], params['J'],
seed=None # important: different disorder each time
)
H = tbrim.build_full_hamiltonian()
for C in C_values:
network = Network()
adj, diag = network.build_network(H, C)
# link counts
link_counts, *_ = network.compute_link_distribution(adj)
# hub selection
hub_index = 0
hub_state = tbrim.basis[hub_index]
occupied_orbitals = tbrim._occupied_orbitals(hub_state)
print(f"Hub state: index={hub_index}, occupied orbitals={occupied_orbitals}")
erdos, reach_frac, *_ = network.compute_erdos_numbers(adj, hub_index)
results[params['label']][C]["all_links"].extend(link_counts)
#results[params['label']][C]["all_erdos"].extend(erdos[erdos >= 0])
results[params['label']][C]["all_erdos"].extend(erdos)
results[params['label']][C]["reach_frac"].append(reach_frac)
# Plot results
figure = PlotFigure()
figure.plot_figure(results, param_sets, C_values)
return results
def analyze_connectivity(self, M, L, C_values):
"""More detailed analysis of network conectivity"""
d = int(comb(M, L, exact=True))
print(f"\nDetailed analysis for M={M}, L={L}, d={d}")
# Create Hamiltionian for two extreme cases
print("\n1. SYK-like (V=0, J=1):")
tbrim_syk = TBRIMHamiltonian(M, L, V=0.0, J=1.0)
H_syk = tbrim_syk.build_full_hamiltonian()
print("\n2. Weak interaction (V=1, J=0.05):")
tbrim_weak = TBRIMHamiltonian(M, L, V=np.sqrt(14), J=0.025)
H_weak = tbrim_weak.build_full_hamiltonian()
# Analyze conectivity
for C in C_values:
print(f"\n{'='*50}")
print(f"Network Analysis with C={C}")
print(f"{'='*50}")
for label, tbrim, H in [("SYK", tbrim_syk, H_syk), ("Weak", tbrim_weak, H_weak)]:
print(f"\n{label}:")
network = Network()
adj, diag = network.build_network(H, C=C)
# Count connections
total_possible = d * (d - 1) // 2
actual_connections = np.sum(adj) // 2
print(f" Total possible connections: {total_possible}")
print(f" Actual connections: {actual_connections}")
print(f" Connection density: {100*actual_connections/total_possible:.2f}%")
# Analyze component structure
components = self.find_connected_components(adj)
print(f" Number of connected components: {len(components)}")
if len(components) > 0:
largest = max(components, key=len)
print(f" Largest component size: {len(largest)}/{d} ({100*len(largest)/d:.1f}%)")
# Compute average degree
degrees = np.sum(adj, axis=1)
print(f" Average degree: {np.mean(degrees):.2f}")
print(f" Max degree: {np.max(degrees)}")
print(f" Min degree: {np.min(degrees)}")
def find_connected_components(self, adj):
"""
Find all connected components in the graph
"""
n = adj.shape[0]
visited = np.zeros(n, dtype=bool)
components = []
for i in range(n):
if not visited[i]:
# BFS to find component
component = []
queue = deque([i])
visited[i] = True
while queue:
node = queue.popleft()
component.append(node)
neighbors = np.where(adj[node])[0]
for neighbor in neighbors:
if not visited[neighbor]:
visited[neighbor] = True
queue.append(neighbor)
components.append(component)
return components
if __name__=="__main__":
# Set random seed for reproducibility
np.random.seed(42)
M = 12 # orbitals
L = 7 # particles
C_values = [0.1, 1.0]
n_realizations=2
# Run
model = RunModelAndAnalysis()
results = model.run_model(M, L, C_values, n_realizations)
# A detailed conectivity analysis
model.analyze_connectivity(M, L, C_values)
The result I get:
which seems prretty much like the Fig9 from the paper (link above). In figure 9 they analyzed the hub=5720 ("corresponding to the many-body state where orbitals 3, 4, 8, 9, 12, 14, 15 out of 16 orbitals are occupied").
I concluded: or I am completely wrong or I am plotting hub=5720 while I was trying to plot hub=0, which means I am messing up with the selection process when I code:
# hub selection
hub_index = 0
hub_state = tbrim.basis[hub_index]
occupied_orbitals = tbrim._occupied_orbitals(hub_state)
print(f"Hub state: index={hub_index}, occupied orbitals={occupied_orbitals}")
erdos, reach_frac, *_ = network.compute_erdos_numbers(adj, hub_index)
If this is the case, that is, I am selecting hub=5720 instead of hub=0, then its a minor problem but I rewrote chunks of this code and the results were worse and worse. This is the best code I have but is not exactly doing what I was trying to do. This is my first time working with something like that, and I am writing this code for several days now.
Any help would be highly appreciated.