{"id":"97db8f73-eda0-46d2-8963-2156c1472e85","arxiv_id":"2505.05361","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"Random boundary illuminations plus a weighted energy estimate yield converged finite element reconstruction rates for quantitative photoacoustic tomography, with errors controlled by noise, mesh size, and regularization.","lead":"This paper proves error estimates for a finite element reconstruction scheme that recovers tissue optical diffusion and absorption coefficients from internal photoacoustic data under random boundary light illuminations. The estimates connect mesh size and regularization strength to data noise, giving practitioners explicit parameter choices.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"No significant objection identified. The central conditional error estimate is coherent; the probabilistic non-zero condition is explicit, and the L2(Ω) claim for D and σ is supported by Theorem 3.1. Only the advertised L-dependence in Remarks 2.4/3.1 is not derived from the stated bounds.","rationale":"The reader's conditional verdict is reasonable. The reader identifies the non-zero condition (2.4) as the weakest assumption; I agree it is the key probabilistic premise, but it is explicitly quantified in Proposition 2.1, so it does not constitute an unacknowledged flaw. The more concrete loose end is the mismatch between the Ω′ bound in Theorem 2.2 and the full-domain L2(Ω) rate claimed in Remarks 2.4 and 3.1. The final reconstruction theorem, Theorem 3.1, does support an L2(Ω) estimate for D and σ because q* is extended outside Ω′ using known coefficients and measured data, with an O(δ) consistency error. However, the specific L^{7/8}δ^{1/4−ε} rate is not a direct consequence of the stated theorems; a direct combination of Proposition 3.1 and Theorem 3.1 yields the same δ-rate with a larger L-power. Since the paper's main scientific claim is the rigorous first-order error estimate in h, η, and δ together with the high-probability framework, this overclaim does not invalidate the central theorem. No code or reproducible artifacts are shipped, which limits independent numerical verification but does not affect the mathematical argument. Overall, the verdict should remain conditional pending correction of the rate statements.","tokens_in":24080,"tokens_out":38691,"duration_ms":409302,"concrete_test":"Verify the final rate by deriving it directly from Proposition 3.1 and Theorem 3.1: set h^2 L^{1/2} ~ δ and α ~ δ^2, and compute η = ||q† − q_h^*||_{L2(Ω′)} from Proposition 3.1. If the sharpest resulting bound for ||D† − D*|| + ||σ† − σ*|| is C L δ^{1/4} rather than C L^{7/8} δ^{1/4−ε}, amend Remarks 2.4 and 3.1 to state the correct L-dependence while keeping the main theorem unchanged.","verdict_should_be":"UNCHANGED","load_bearing_attack":"After checking the chain from Proposition 2.1 through Theorem 2.2 to Theorem 3.1, I do not find a load-bearing flaw in the central claim. The non-zero condition (2.4) is genuinely the enabling premise, but it is stated as a high-probability event with explicit probability (2.3), so the theorem's conditional nature is not a hidden assumption. The intermediate IDP estimate is only on Ω′, but Theorem 3.1 handles the outside region by setting q* = D†(Zδ^(1)/σ†)^2 there, yielding an O(δ) contribution, so the abstract's L2(Ω) claim for D and σ is supported. The only soft spot is in Remarks 2.4 and 3.1: the advertised L^{7/8}δ^{1/4−ε} rate is not derived from Theorem 2.2 and Proposition 3.1; the bound that follows from those results is O(L δ^{1/4}) with the same δ-exponent. This affects the stated L-dependence, not the existence of the convergence estimate.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper develops and analyzes a two-stage finite element method for quantitative photoacoustic tomography (QPAT) in the diffusive regime, reconstructing the diffusion coefficient D and absorption coefficient σ from internal deposited-energy data generated by random boundary illuminations. In the first stage, the problem is reduced to an inverse diffusivity problem (IDP) for q = D u_1^2, which is solved by a regularized output least-squares formulation with P1 finite elements. In the second stage, a direct elliptic problem for v = 1/u_1 - 1 is solved and D, σ are recovered algebraically. The main theoretical results are: a high-probability non-zero gradient condition under random boundary data (Proposition 2.1), a conditional stability estimate (Theorem 2.1), an L2(Ω′) error estimate for the discrete diffusivity (Theorem 2.2), and an L2(Ω) error estimate of order h+η+δ for the final coefficients (Theorem 3.1), with the parameter choice h^2 L^{1/2} ∼ δ and α ∼ δ^2 leading to a δ^{1/4} convergence rate. Numerical experiments on smooth and nonsmooth coefficients illustrate the predicted behavior.","tokens_in":24326,"tokens_out":17077,"duration_ms":161465,"significance":"If the estimates are correct, this is a substantial contribution: it provides a rigorous finite element error analysis for QPAT with randomly chosen illuminations, a regime in which most existing analyses are at the continuous level. The proof chain is coherent and uses appropriate tools: the published probabilistic non-zero condition of [1], the weighted energy estimate of [13,30], and the decoupling of QPAT into an inverse diffusivity problem followed by a direct solve. The authors are explicit about the probabilistic nature of the non-zero condition and about the fact that the IDP estimate is on Ω′. The paper also gives concrete guidance for selecting the mesh size and regularization parameter from the noise level, and the numerical rates are consistent with the predicted δ^{1/4} behavior. These strengths make the paper valuable to the numerical analysis and inverse problems communities.","major_comments":[{"comment":"The L2(Ω) convergence rate advertised in Remarks 2.4 and 3.1 is not a direct consequence of the stated theorems. Theorem 2.2 proves an L2(Ω′) estimate for the diffusivity, and substituting h^2 L^{1/2} ∼ δ and α ∼ δ^2 into that estimate gives ∥q†−q*_h∥_{L2(Ω′)} ≤ C L^{7/8}δ^{1/4} for d=2 and C L^{(7+ε)/8}δ^{(1−ε)/4} for d=3, not the stated C L^{7/8}δ^{1/4−ε}. Moreover, the interpolating argument in Remark 2.4 uses H2 regularity for w(q) with q replaced by the discrete reconstruction q*_h, with constants independent of h; this regularity is not established for piecewise-linear q*_h. Since Remark 3.1 repeats the same rate for the final D and σ, the claim should be either proved from Theorem 2.2 and Theorem 3.1 with correct exponents, or explicitly labeled as heuristic.","section":"§2.2 (Theorem 2.2, Remark 2.4) and §3 (Remark 3.1)"},{"comment":"The abstract states that the paper provides 'a rigorous error estimate in L2(Ω) norm for the numerical reconstruction' without specifying that the intermediate diffusivity estimate is proved only on Ω′. This distinction matters because the IDP solver itself is not shown to be accurate in L2(Ω); the full-domain result is obtained in Theorem 3.1 only after setting q* outside Ω′ to the known coefficient value. Please state the Ω′/Ω distinction explicitly in the abstract and in Remark 2.4 to avoid overclaiming the scope of the IDP result.","section":"Abstract and §2.2"}],"minor_comments":[{"comment":"The manuscript states 'Ω ⊂ R^d (d = 2, 2)'; this should read '(d = 2, 3)'.","section":"§2.2, first paragraph"},{"comment":"The text refers to 'Assumption 3.1(iii)', but Assumption 3.1 has only items (i) and (ii); the condition g^(1) ≡ 1 appears in item (ii). Please correct the cross-reference.","section":"§4.1"},{"comment":"In several displayed estimates the proof writes 'u(ℓ)' where the intended quantity is 'w(ℓ)(q†)'; please make the notation consistent.","section":"Lemma 2.2 proof"},{"comment":"For nonhomogeneous Dirichlet problems on polygonal domains, H2 regularity requires compatibility conditions on the boundary data. Please state the required condition or cite the precise theorem that covers the present case.","section":"Remark 2.5"},{"comment":"The numerical experiments do not describe the optimization solver used for the least-squares problem, the initialization, or the random seeds for the boundary illuminations. Reporting these details would improve reproducibility; the reported convergence rates also appear to come from single realizations.","section":"§4"},{"comment":"The piecewise-constant coefficients in Example 4.5 do not satisfy Assumption 3.1; the favorable reconstructions there should be framed as additional numerical evidence rather than as validation of the theory.","section":"Example 4.5"}],"recommendation":"major_revision","confidential_remarks":"The core proof chain appears sound, and the conditional estimates are clearly stated. The main issue is an overstatement in the rate claims in Remarks 2.4 and 3.1, which should be corrected or downgraded. With those changes, this would be a solid contribution to the numerical analysis of QPAT."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: this is a solid numerical analysis paper and the main theorem is probably right. It is the first FEM error analysis I know for the two-parameter QPAT problem in a diffusive regime with a vanishing source and product measurement H = σu. The trick is to combine the random boundary illuminations of [1] with weighted energy estimates from [30, 18] in a two-stage least-squares scheme: first solve an inverse diffusivity problem for q = D u_1^2, then solve a direct problem for 1/u_1 and recover D and σ. Theorem 3.1 gives ∥D† - D*∥_{L2(Ω)} + ∥σ† - σ*∥_{L2(Ω)} ≤ C(h + η + δ) with η the diffusivity error, and Proposition 3.1 gives the diffusivity error on Ω′. I checked the chain from Proposition 2.1 through Theorem 2.2 to Theorem 3.1; I do not find a load-bearing flaw. The non-zero condition is genuinely probabilistic, but it is explicit with probability (2.3), so the conditional nature is not hidden.\n\nThe numerical section is more than decorative: four smooth/nonsmooth examples plus a piecewise-constant case outside the theory, and observed rates align with the predicted δ^{1/4} (some slightly better). No code or data are shipped, so reproducibility is limited but not a fatal issue for this type of paper.\n\nSoft spots, in decreasing order. (1) The abstract and Section 2 claim an L2(Ω) error estimate for the diffusivity reconstruction, but Theorem 2.2 proves the bound only on Ω′. The final QPAT theorem recovers an L2(Ω) statement for D and σ by defining q* outside Ω′ using the datum Zδ^(1)/σ†, which gives an O(δ) contribution; so the abstract's claim is supported at the end, but the intermediate L2(Ω) diffusivity claim is not. (2) The advertised L^{7/8}δ^{1/4-ε} rate in Remarks 2.4 and 3.1 is not actually derived from the stated bounds; what follows is O(L δ^{1/4}) with the same δ exponent, so the L-dependence is weaker than advertised. This matters for users who want to choose L, but it does not break the convergence result. (3) The proof of Theorem 2.2 relies on elliptic regularity and L2-projection stability in Lq for q < 2d/(d-2); the estimates are sketched but I did not find an inconsistency.\n\nWho is this for? Researchers in quantitative photoacoustic tomography and inverse coefficient identification, especially those using FEM or random-illumination strategies. It deserves a serious referee. I would recommend acceptance after the authors correct the domain mismatch and either prove or demote the advertised L-dependence.","headline":"Solid first FEM error analysis for two-parameter QPAT with vanishing source; the main theorem holds, but the advertised L2(Ω) and L-dependence overstate what is proved.","tokens_in":24843,"tokens_out":2045,"would_cite":true,"duration_ms":19105,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65N21","65N30","35R30"],"pacs":[],"model":"deepseek-v4-flash","headline":"A two-stage finite element scheme recovers the diffusion and absorption coefficients in quantitative photoacoustic tomography with L2 error of order h+η+δ, on a high-probability non-zero-gradient event.","keywords":["quantitative photoacoustic tomography","inverse diffusivity problem","random boundary illuminations","finite element approximation","weighted energy estimate","least-squares regularization","convergence rate"],"falsifier":"On a fixed smooth coefficient pair, compute the quantity $\\max_{\\ell}|\\nabla w^{(\\ell)}\\cdot\\nu|$ on a fine grid over $\\Omega'$ for many random illumination draws; on the draws where it stays below the $C_0/2$ threshold on a positive-measure set, run the two-stage scheme with noise-free data and check whether the $L^2$ error still decays at the predicted rate in $h$: if it does, the non-zero gradient condition is not necessary, and if it does not, the theorem's key premise is confirmed.","tokens_in":23862,"feed_emoji":"🩺","tokens_out":12467,"duration_ms":114819,"temperature":0.7,"pith_summary":"The paper claims that quantitative photoacoustic tomography in the diffusive regime can be turned into a provably convergent numerical scheme: recover the diffusivity coefficient $q = D u_1^2$ from ratios of internal energy measurements via a regularized least-squares finite element method, then recover $u_1$ from a direct elliptic solve and read off $D$ and $\\sigma$ algebraically. The central result is a rigorous $L^2(\\Omega)$ error bound for both coefficients of order $h + \\eta + \\delta$, where $h$ is the mesh size, $\\delta$ the data noise, and $\\eta$ the error of the first-stage diffusivity reconstruction. The bound holds with high probability, because random boundary illuminations make the gradient of the quotient solutions non-zero on the reconstruction region. A reader should care because the analysis gives an explicit rule for choosing mesh size and regularization parameter from the measured noise level, turning a heuristic two-step pipeline into a scheme with a proven rate.","feed_headline":"Photoacoustic coefficient recovery gets first-order error bounds","feed_subtitle":"Random boundary lights stabilize the inversion and guide mesh and regularization choices from noise.","key_machinery":"The load-bearing object is the quotient $w^{(\\ell)}=H^{(\\ell+1)}/H^{(1)}=u^{(\\ell+1)}/u^{(1)}$, which satisfies the one-parameter elliptic equation $-\\nabla\\cdot(q\\nabla w)=0$ with $q=D u_1^2$. The main estimate is a weighted energy identity obtained by testing the weak form with $\\varphi=(q-\\tilde q)w/q$; it puts $\\frac12\\int_{\\Omega'}\\frac{(q-\\tilde q)^2}{q}|\\nabla w|^2\\,dx$ on the left-hand side, so the non-zero gradient condition (2.4) turns a small misfit in the quotient data into a small $L^2$ error in $q$. The non-zero condition itself is supplied probabilistically: boundary illuminations are drawn as Gaussian series in an $H^{1/2}(\\partial\\Omega)$ orthonormal basis, and Proposition 2.1 shows that with probability at least $1-L^d e^{-C_1 L}-L e^{-C_2 M}$ some directional derivative of the quotient solutions is bounded away from zero on $\\Omega'$. The numerical side is a standard piecewise-linear Galerkin method with an $H^1$-seminorm penalty, and the discrete error analysis combines interpolation, inverse, and duality estimates to transfer the continuous stability bound to the finite element solution.","core_discovery":"The central discovery, stated as Theorem 3.1, is that the two-stage procedure reconstructs both optical coefficients with $L^2(\\Omega)$ error of order $h+\\eta+\\delta$. In the first stage the quotient $w^{(\\ell)} = H^{(\\ell+1)}/H^{(1)}$ is an observed function satisfying $-\\nabla\\cdot(q^\\dagger \\nabla w^{(\\ell)})=0$, so the paper recovers $q^\\dagger = D^\\dagger|u^{(1)}|^2$ by minimizing a regularized least-squares misfit over piecewise-linear finite elements; Theorem 2.2 controls this first-stage error, and balancing $h^2 L^{1/2}\\sim\\delta$ with $\\alpha\\sim\\delta^2$ yields the rate $L^{7/8}\\delta^{1/4-\\epsilon}$ in two dimensions. In the second stage, $v=1/u^{(1)}-1$ solves the direct problem $-\\nabla\\cdot(q^\\dagger\\nabla v)=H^{(1)}$ with zero boundary data, and replacing $q^\\dagger$ and $H^{(1)}$ by their numerical counterparts gives $v_h$, from which $D^* = q^*|v_h+1|^2$ and $\\sigma^* = Z_\\delta^{(1)}(v_h+1)$ are formed. The proof uses a weighted energy identity with the special test function $\\varphi=(q-\\tilde q)w/q$, which converts the non-zero gradient condition into a lower bound on the data misfit and hence into an upper bound on the coefficient error. All statements hold with the probability in (2.3), and the error constant is independent of $h$, $\\delta$, and $\\alpha$.","pith_inferences":["Since $w^{(\\ell)}$ is formed directly from the measured energies, the non-zero gradient condition could be monitored on a computational grid before inverting, enabling an adaptive acquisition rule that adds random illuminations until condition (2.4) is observed rather than relying only on the probabilistic guarantee.","The $\\delta^{1/4}$ rate is probably an artifact of the $L^2$ misfit and $H^1$ penalty; using the weighted energy structure itself as the data fidelity term could restore the $\\delta^{1/2}$ rate suggested by the conditional stability estimate.","The same quotient reduction to a source-free inverse diffusivity problem should transfer to other hybrid imaging modalities with boundary-only illumination and internal data, such as conductivity or fluorescence imaging.","The experiments with piecewise-constant coefficients suggest the smoothness assumptions are sufficient but not necessary; testing on L-shaped domains or discontinuous coefficients would map the true boundary of validity."],"forward_implications":["Choosing $h^2 L^{1/2}\\approx\\delta$ and $\\alpha\\approx\\delta^2$ gives a predicted rate of order $\\delta^{1/4-\\epsilon}$ for both coefficients in two dimensions, and the numerical experiments report exponents between $0.22$ and $0.42$ for the relative $L^2$ errors.","The theorem supplies a parameter-selection rule: mesh size and regularization strength should be scaled as $h\\sim\\delta^{1/2}$ and $\\alpha\\sim\\delta^2$ once the noise level is known.","The scheme covers the practical case of vanishing source and boundary-only illumination, which earlier two-observation reconstructions could not handle because their positivity conditions failed.","Because the final error is linear in the first-stage diffusivity error $\\eta$, any improvement in numerical inverse diffusivity solvers transfers directly to quantitative photoacoustic tomography.","The piecewise-constant and non-smooth experiments still converge at the predicted rates, indicating that the stated regularity assumptions are sufficient but likely not necessary."],"supporting_citations":[{"why":"Supplies the probabilistic non-zero gradient condition used as Proposition 2.1, the key stability input for all subsequent bounds.","marker":"[1]"},{"why":"Introduces the decoupling quotient w=u2/u1=H2/H1 that reduces QPAT to an inverse diffusivity problem.","marker":"[12]"},{"why":"Extends the decoupled procedure to the multi-source diffusive regime and provides uniqueness and Hölder stability for QPAT.","marker":"[11]"},{"why":"Develops the weighted energy estimate with a special test function that the paper adapts to the vanishing-source, boundary-illumination setting.","marker":"[30]"},{"why":"Provides the weighted-energy stability framework that motivates the discrete error analysis.","marker":"[13]"},{"why":"Gives the earliest discrete least-squares scheme for the inverse diffusivity problem and the error-estimate template O(h^r+h^{-2}delta).","marker":"[22]"},{"why":"Supplies improved finite element error estimates for parameter identification under stronger non-zero conditions.","marker":"[43]"},{"why":"Recent two-parameter reconstruction from two internal measurements whose two-step structure the present method extends to random illuminations.","marker":"[18]"}],"fun_headline_variants":["First-order error bounds for photoacoustic coefficient recovery","Random boundary lights stabilize photoacoustic inversion","Finite element method yields guaranteed error for photoacoustic tomography","Reconstructing optical coefficients with provable accuracy","Stable photoacoustic inversion via random boundary illumination"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"All the error bounds collapse if the random boundary illuminations do not make the quotient solutions satisfy $\\max_{\\ell}|\\nabla w^{(\\ell)}\\cdot\\nu|\\ge C_0$ on $\\Omega'$, and the paper proves this only with overwhelming probability, not deterministically.","fun_headline_variants_meta":{"raw":{"variants":["First-order error bounds for photoacoustic coefficient recovery","Random boundary lights stabilize photoacoustic inversion","Finite element method yields guaranteed error for photoacoustic tomography","Reconstructing optical coefficients with provable accuracy","Stable photoacoustic inversion via random boundary illumination"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.00108,"raw_usage":{"total_tokens":4585,"prompt_tokens":1078,"completion_tokens":3507,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":694,"completion_tokens_details":{"reasoning_tokens":3432}},"tokens_in":694,"tokens_out":3507,"duration_ms":24489,"temperature":1.0,"reasoning_tokens":3432,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T23:06:04.297792+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a fixed smooth coefficient pair, compute the quantity $\\max_{\\ell}|\\nabla w^{(\\ell)}\\cdot\\nu|$ on a fine grid over $\\Omega'$ for many random illumination draws; on the draws where it stays below the $C_0/2$ threshold on a positive-measure set, run the two-stage scheme with noise-free data and check whether the $L^2$ error still decays at the predicted rate in $h$: if it does, the non-zero gradient condition is not necessary, and if it does not, the theorem's key premise is confirmed.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the probabilistic non-zero gradient condition used as Proposition 2.1, the key stability input for all subsequent bounds."},{"cited_title":"Bal and G","cited_arxiv_id":null,"evidence_quote":"Introduces the decoupling quotient w=u2/u1=H2/H1 that reduces QPAT to an inverse diffusivity problem."},{"cited_title":"Bal and K","cited_arxiv_id":null,"evidence_quote":"Extends the decoupled procedure to the multi-source diffusive regime and provides uniqueness and Hölder stability for QPAT."},{"cited_title":"Jin and Z","cited_arxiv_id":null,"evidence_quote":"Develops the weighted energy estimate with a special test function that the paper adapts to the vanishing-source, boundary-illumination setting."},{"cited_title":"Bonito, A","cited_arxiv_id":null,"evidence_quote":"Provides the weighted-energy stability framework that motivates the discrete error analysis."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Gives the earliest discrete least-squares scheme for the inverse diffusivity problem and the error-estimate template O(h^r+h^{-2}delta)."},{"cited_title":"Wang and J","cited_arxiv_id":null,"evidence_quote":"Supplies improved finite element error estimates for parameter identification under stronger non-zero conditions."},{"cited_title":"Cen and Z","cited_arxiv_id":null,"evidence_quote":"Recent two-parameter reconstruction from two internal measurements whose two-step structure the present method extends to random illuminations."}],"review_version":1}