コンテンツにスキップ

2025年度

25/10/04 進捗報告

SNonlyであってもCviscパラメータは関係あるので、すでにあるSNonlyだと200πになってしまっている。

そこで、100πに変更してrunを回すことにした。

[z6b481@squidhpc3 L50N512_Cvisc100piSNonly]$ qstat
RequestID       ReqName  UserName Queue     Pri STT S   Memory      CPU   Elapse R H M Jobs
--------------- -------- -------- -------- ---- --- - -------- -------- -------- - - - ----
897746.sqd      L50N512C z6b481   SC64        0 QUE -    0.00B     0.00        0 Y Y Y   64

ToDo-1.3

Χ^2を計算する

計算条件

  • 比較対象: McConnell & Ma (2013) の \(z \approx 0\) の BH–Mstar 関係
  • 使用データ: 観測点 32 個
  • 誤差: 観測誤差（σ）を使用、シミュレーションの分位幅は使用せず
  • 評価指標: reduced χ² = χ² / N

結果

モデル名 reduced χ² 備考
\(C_\text{visc} = 100\pi\) 8.19 観測に最も近い
\(C_\text{visc} = 200\pi\) 16.53 ややずれが大きい
SNonly 53.86 最も大きく乖離

※このSNonlyは、\(C_{\rm visc} = 200\pi\) である。

BH mass–Star mass / BH mass–Stellar surface density in 3 kpc

image-20251008212222106

変化点を真面目に推定する

1 変化点の折れ線回帰（hinge/broken-stick）で “ぽきっ” の位置を推定し、ブートストラップで信頼区間も出した。

使用したコード

# Fit a 1-break "hinge" (broken-stick) regression y = c + m*x + k*max(0, x-xb)
# using grid search over xb + closed-form least squares for (c,m,k).
import numpy as np, pandas as pd, json, os
import matplotlib.pyplot as plt

df = pd.read_csv("/mnt/data/xx.csv")
x = df["x"].to_numpy(dtype=float)
y = df["y"].to_numpy(dtype=float)

def fit_given_xb(xb, x, y, w=None):
    X = np.column_stack([
        np.ones_like(x),
        x,
        np.maximum(0.0, x - xb),
    ])
    if w is None:
        beta, *_ = np.linalg.lstsq(X, y, rcond=None)
        resid = y - X @ beta
    else:
        W = np.diag(w)
        beta = np.linalg.inv(X.T @ W @ X) @ (X.T @ W @ y)
        resid = y - X @ beta
    rss = float((resid**2).sum())
    return beta, rss

# Coarse-to-fine search for xb
xmin, xmax = float(np.min(x)), float(np.max(x))
grid1 = np.linspace(xmin, xmax, 401)
rss1 = []
for xb in grid1:
    _, rss = fit_given_xb(xb, x, y)
    rss1.append(rss)
xb1 = grid1[int(np.argmin(rss1))]

# refine around xb1
span = (xmax - xmin) * 0.1 if (xmax - xmin) > 0 else 0.1
grid2 = np.linspace(max(xmin, xb1 - span), min(xmax, xb1 + span), 401)
rss2 = []
for xb in grid2:
    _, rss = fit_given_xb(xb, x, y)
    rss2.append(rss)
xb_best = grid2[int(np.argmin(rss2))]
beta_best, rss_best = fit_given_xb(xb_best, x, y)
c, m, k = beta_best.tolist()

# Bootstrap CI for xb
rng = np.random.default_rng(0)
B = 1000
xb_samples = []
n = len(x)
for _ in range(B):
    idx = rng.integers(0, n, n)
    xb1s = np.linspace(xmin, xmax, 201)
    # coarse
    best = None; best_rss = np.inf
    for xb in xb1s:
        _, rss = fit_given_xb(xb, x[idx], y[idx])
        if rss < best_rss:
            best_rss = rss; best = xb
    # fine
    span = (xmax - xmin) * 0.1 if (xmax - xmin) > 0 else 0.1
    xb2s = np.linspace(max(xmin, best - span), min(xmax, best + span), 201)
    best2 = None; best_rss2 = np.inf
    for xb in xb2s:
        _, rss = fit_given_xb(xb, x[idx], y[idx])
        if rss < best_rss2:
            best_rss2 = rss; best2 = xb
    xb_samples.append(best2)

xb_samples = np.array(xb_samples)
xb_ci = (float(np.percentile(xb_samples, 2.5)), float(np.percentile(xb_samples, 97.5)))
xb_std = float(np.std(xb_samples, ddof=1))

# Save results
res = {
    "break_xb": float(xb_best),
    "break_ci_95_lower": xb_ci[0],
    "break_ci_95_upper": xb_ci[1],
    "break_std_bootstrap": xb_std,
    "intercept_c": float(c),
    "slope_left_m": float(m),
    "slope_right_m_plus_k": float(m + k),
    "k_delta_slope": float(k),
    "rss": float(rss_best),
    "n": int(n),
}

res_df = pd.DataFrame([res])
res_path = "/mnt/data/break_fit_results.csv"
res_df.to_csv(res_path, index=False)

# Plot: data + fit + break with CI band
xs = np.linspace(xmin, xmax, 400)
ys = c + m*xs + k*np.maximum(0.0, xs - xb_best)

plt.figure(figsize=(6,4))
plt.scatter(x, y, s=10, alpha=0.5)
plt.plot(xs, ys, linewidth=2)
plt.axvline(xb_best, linestyle="--")
plt.fill_betweenx([min(y), max(y)], xb_ci[0], xb_ci[1], alpha=0.15)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Broken-stick fit with estimated breakpoint")
plt.tight_layout()
plot_path = "/mnt/data/broken_stick_fit.png"
plt.savefig(plot_path, dpi=160)
plt.close()

res_path, plot_path, res

推定結果(left panel)

  • 変化点 \(x_b = 9.54^{+0.12}_{-0.18}\)
  • 左側の傾き \(m = 0.079\)
  • 右側の傾き \(m+k = 2.197\)
  • 切片 \(c =4.559\)
  • 図（データ点＋フィット線、縦の破線が変化点、半透明帯が95%CI）

broken_stick_fit

推定結果(right panel)

  • 変化点 \(x_b = 8.09^{+0.43}_{-0.18}\)
  • 左側の傾き \(m = 0.071\)
  • 右側の傾き \(m + k = 2.436\)
  • 傾き差 \(k = 2.366\)
  • 切片 \(c = 4.737\)
  • サンプル数 \(n = 32\)
  • 図（データ点＋フィット線＋変化点＋95%信頼区間）

broken_stick_fit_xx2

mention

Recent cosmological zoom-in simulations, such as the FIRE project \citep{byrneStellarFeedbackregulatedBlack2023}, have revealed that black hole (BH) growth proceeds through two distinct phases.  
In the early phase, stellar feedback from supernova explosions and stellar winds repeatedly expels gas from the galactic center.  
This process prevents a stable inflow of material onto the black hole, resulting in weak and intermittent accretion.  
As the galaxy evolves and its gravitational potential deepens, the situation changes substantially.  
The galactic disk settles into a stable and thin configuration, while the surrounding gas becomes thermally supported.  
This thermal support suppresses further gas outflows driven by feedback.  
As a result, a long-lived gas reservoir forms in the central region, leading to a transition toward a sustained and rapid BH growth phase.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=\linewidth]{pic/fig_bh_stellar_apj.pdf}
    \caption{Caption}
    \label{fig:bh-stellar-dis}
\end{figure}

We observe a similar two-phase behavior in our simulations.  
When examining the median BH growth histories across different halo mass bins, we find that this transition occurs almost independently of halo mass.  
Figure \ref{fig:bh-stellar-dis} shows these relations: the left panel presents BH mass versus stellar mass, and the right panel shows BH mass versus stellar surface density within 3 kpc.  
For each redshift, we compute the median and the 16th–84th percentile range of BH and stellar mass, and plot them with error bars.  
The solid lines connect the median values and represent the statistical evolutionary tracks.  
In both panels, regardless of halo mass, a clear transition point can be identified where the growth mode changes.

To quantify this break, we applied a single-break (broken-stick) regression model of the form  $y = c + mx + k\max(0,x - x_b)$,  
where $x_b$ represents the breakpoint, $m$ is the slope below the break, and $m + k$ is the slope above it.  
We determined $x_b$ by performing a grid search that minimizes the sum of squared residuals, and estimated the other parameters $(c, m, k)$ using linear least squares.  
The uncertainty of the breakpoint was evaluated by 1000 bootstrap resamplings.  
The best-fit breakpoint was found to be $x_b = 9.54^{+0.12}_{-0.18}$, indicating a distinct change in slope around $\log_{10} M_{\mathrm{stellar}} \simeq 9.54$.