{"id":"26423afa-4df8-4894-b198-5d8d65de6041","arxiv_id":"2608.03879","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"high","formal_verification":"none","parameter_count":5,"one_line_summary":"Closed-form expressions for arbitrary-order multi-time correlation functions with Duschinsky rotation and finite temperature are derived and applied to 2D resonance Raman spectra.","lead":"This paper derives closed-form equations for computing high-order nonlinear spectroscopy signals, including molecular vibrations that rotate and change frequency between electronic states, at finite temperature. It applies them to simulate fifth-order two-dimensional resonance Raman spectra, showing that such spectra carry fingerprints of vibrational rotations.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (27) for the linear coefficient E is inconsistent with the Gaussian integral (Eq. 16) for n≥3 and with the multi-state formula (Eq. 36), so the central two-state response expression (Eq. 29) is unreliable as printed.","rationale":"The reader identified the preresonance/resonance mismatch as the weakest assumption and mentioned Eq. (27) only in passing. In my reading, the Eq. (27) inconsistency is the more load-bearing concern because it challenges the core closed-form derivation (Eq. 29) that underlies the entire framework. If the code follows Eq. (27), all 2DRR spectra are wrong; if it follows Eq. (36), the manuscript still contains a central equation that is incorrect as printed and must be corrected or clarified. The preresonance mismatch is a real application-level problem, but it would only affect the illustrative spectra, not the general method's validity. Since both issues require correction before the results can be fully trusted, the reader's CONDITIONAL verdict is unchanged, though for a different primary reason. The proposed test (direct numerical comparison or code modification) would settle whether the Eq. (27) error is merely typographical or actually affects the published spectra.","tokens_in":18233,"tokens_out":21245,"duration_ms":178504,"concrete_test":"For N=1, J=1, K=1, take n=5 and fixed τ0...τ5 (e.g., all equal to 1). Compute C(τ0,...,τ5) two ways: (1) numerically integrate Eq. (16) directly (or evaluate the Gaussian integral with the correct E derived from Eq. (16)); (2) evaluate Eq. (29) using Eq. (27) for E. If the results differ, Eq. (27) is invalid. Additionally, run the published GitHub code for model 3 in Table 2 with Eq. (27) and with the corrected E from Eq. (36) specialized to two states; if the 2DRR spectra change, the paper's central results depend on an incorrect printed expression.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The two-state derivation's linear coefficient E, Eq. (27), does not follow from the Gaussian integral in Eq. (16) for any n>1. For odd n (as in the fifth-order 2DRR application), the ground-state propagators occur at even time indices τ0, τ2, ..., τ_{n−1}. Expanding Eq. (16) into (i/2)z^T D z + i E^T z assigns the displacement K to coordinates via the ground-state Gaussian factors: the propagator at τ_{2k} (k≥1) couples z_{2k−1} and z_{2k}, giving E_{2k−1}=E_{2k}=J^T c(τ_{2k})K; the propagator at τ0 couples z0 and z_n, giving E0=En=J^T c(τ0)K. Equation (27) instead pairs E_i with c(τ_{⌊i/2⌋}), which for n=3 gives (c0,c0,c1,c1) instead of the correct (c0,c2,c2,c0); for n=5 it gives (c0,c0,c1,c1,c2,c2) instead of (c0,c2,c2,c4,c4,c0). The multi-state expression, Eq. (36), reduces to the correct form when specialized to two states. Thus Eq. (29) is wrong as printed; the numerical 2DRR spectra are trustworthy only if the provided code implements the multi-state formula rather than the printed two-state one. This is the most load-bearing concern because it directly undermines the paper's central derivation, not merely an application-level approximation.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript derives closed-form expressions for multi-time correlation functions in harmonic systems with Duschinsky rotation, frequency changes, displacements, and finite temperature, with polynomial scaling in the number of modes. The central formal result is Eq. (29) for a two-state alternating sequence, with a general multi-state extension in Eq. (33). The authors apply the method to fifth-order two-dimensional resonance Raman (2DRR) spectra of two-mode model systems and naphthalene, comparing displaced-harmonic, curvature-change, and full-Duschinsky models. They report that Duschinsky rotation changes the relative intensities of diagonal and cross peaks, increases cross-peak asymmetry, and enhances difference-frequency features, concluding that 2DRR is particularly sensitive to normal-mode rotation between electronic states.","tokens_in":18702,"tokens_out":17474,"duration_ms":153049,"significance":"If correct, the method is a valuable contribution: it extends exact harmonic response theory to arbitrary order without an explicit sum over vibrational eigenstates, handles Duschinsky rotation and finite temperature, and is accompanied by open-source code and ab initio input data. The multi-state formula Eq. (33) appears internally consistent, and the numerical spectra provide concrete, testable predictions. However, the two-state presentation contains a concrete algebraic error in Eq. (27), and the application-level pathway reduction is justified by a preresonance assumption while the calculations are performed at resonance. Because these issues are central to the derivation as printed and to the spectroscopic interpretation, the manuscript requires revision before the results can be fully relied upon.","major_comments":[{"comment":"Equation (27) does not give the correct linear coefficient E for the Gaussian integral (16) when n>1. Expanding the ground-state factors G(Jq+K,Jq'+K,a_g,b_g) assigns a linear term J^T c(τ_i)K to both endpoints of each ground-state propagator. For odd n (e.g., n=3), the ground-state propagators occur at even times τ0,τ2,...; the correct E-vector is (J^T c(τ0)K, J^T c(τ2)K, J^T c(τ2)K, J^T c(τ0)K), whereas Eq. (27) gives (J^T c(τ0)K, J^T c(τ0)K, J^T c(τ1)K, J^T c(τ1)K). The multi-state formula (36) reduces to the correct form when specialized to two states, so Eq. (29) is inconsistent with Eq. (36) as printed. This is a load-bearing error in the central two-state derivation and must be corrected; the code should be checked against the corrected expression.","section":"Sec. 2.1, Eq. (27)"},{"comment":"The 2DRR calculation retains only four of the sixteen fifth-order Liouville pathways, justified by the statement that 'these pathways could be selected using the preresonance regime.' However, Sec. 2.5 states that all spectra are computed with the carrier frequency set to the adiabatic excitation energy, ω_L=ω_eg, i.e., under resonant excitation. Under resonant conditions the omitted twelve pathways are not generically negligible. Unless the authors demonstrate that the selected four pathways dominate at resonance, or recompute with the full set of pathways, the predicted Duschinsky signatures—cross-peak asymmetry, difference peaks, relative intensity changes—may not correspond to the actual 2DRR signal. The excitation regime should be stated explicitly and the pathway reduction justified for the parameters actually used.","section":"Sec. 2.4 and Sec. 2.5"},{"comment":"The two-state derivation is written as if it holds for arbitrary order n, but Eq. (5) is only consistent when n is odd: the sequence of electronic Hamiltonians alternates H_g, H_e, H_g, H_e, ... and the final operator is H_e only for odd n. For even n the last propagator would be H_g, and the D-matrix definitions (20)–(21) would not apply. The applications here use n=1 and n=5, so this does not affect the numerical results, but the 'arbitrary order' claim in the abstract should be qualified, or the two-state formulas should be generalized to arbitrary n.","section":"Sec. 2.1, Eq. (5)"}],"minor_comments":[{"comment":"The pulse duration Δt is not specified numerically, although the approximations t2≈T1 and t4≈T2 rely on Δt being short compared to the vibrational periods. Please state the values used and confirm the validity of the short-pulse assumption for the presented spectra.","section":"Sec. 2.5"},{"comment":"The figure caption lists panels c and d for the antisymmetric components, but the text refers to Fig. 6d for the full-model antisymmetric component. Please verify the panel labels and their correspondence in the text.","section":"Figure 6"},{"comment":"The title contains a spacing artifact ('T emperature'); also, the abstract's 'arbitrary-order' claim should be qualified in view of the odd-n restriction in the two-state derivation, or the general multi-state formula should be highlighted as the arbitrary-order result.","section":"Abstract and title"}],"recommendation":"major_revision","confidential_remarks":"The Eq. (27) inconsistency is the main technical concern; it is a printed error in the two-state derivation. Because Eq. (36) is correct and the code is publicly available, the issue is likely fixable in revision. I would also encourage the authors to verify the pathway reduction by a one-off full-pathway calculation at resonance, as the preresonance justification is not matched by the computational setup."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know. First, this is the first practical closed-form formalism I've seen that handles arbitrary-order vibronic response functions with Duschinsky rotation and finite temperature without a sum over states, and the authors back it with code and data. Second, the printed two-state formula for the linear coefficient E (Eq. 27) is wrong for n>1, and the 2DRR simulations are run at exact resonance while the pathway-selection argument assumes preresonance; both need to be addressed before the central spectra are fully trustworthy.\n\nThe path-integral derivation itself is a genuinely nice piece of work. Reducing the multi-time correlation function to a single Gaussian integral and then solving it by matrix operations is elegant, and the general multi-state expression (Eq. 33) is internally consistent. The application to 2DRR, with the model systems and naphthalene, produces a concrete physical message: Duschinsky rotation shows up mainly as an asymmetry in cross-peak intensities rather than in peak positions. That is useful and well-illustrated. The code is on GitHub and the data on Zenodo, which is how it should be.\n\nThe stress-test about Eq. (27) lands. The correct pattern, as the multi-state Eq. (36) confirms, is E=(c0,c2,c2,c0) for n=3, not (c0,c0,c1,c1). I suspect the code uses the multi-state formula, so the numerical results may well be fine, but as printed the two-state derivation is inconsistent with the Gaussian integral. That is a load-bearing typo, not a cosmetic one.\n\nThe resonance/preresonance point is also real. The paper justifies keeping only 4 of 16 pathways by invoking the preresonance regime, then sets ω_L to the adiabatic excitation energy for all 2DRR calculations. At exact resonance, other Liouville pathways are not obviously negligible. The authors should either run at a detuned frequency or argue why the missing pathways drop out at resonance.\n\nMinor: the parity convention in Eq. (5) is ambiguous for even n, but the interesting cases (n=3,5) are odd, so I'd only ask for a footnote.\n\nWho benefits: anyone computing or interpreting high-order vibrational spectroscopies, especially 2DRR and related techniques. The formalism is clearly presented, and the code is a practical asset.\n\nMy recommendation: send it to review, but ask for a corrected Eq. (27), an explicit statement of which E the code uses, and a thoughtful response on the preresonance assumption. If those are fixed, it's a solid contribution.","headline":"Solid arbitrary-order harmonic response-function formalism with one printed equation error and a resonance-vs-preresonance mismatch that need fixing before the 2DRR results are fully persuasive.","tokens_in":19164,"tokens_out":7001,"would_cite":true,"duration_ms":56465,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper derives closed-form, polynomial-scaling equations for arbitrary-order vibronic response functions with Duschinsky rotation and finite temperature, and uses them to identify distinct Duschinsky signatures in 2D resonance Raman spe","keywords":["Multi-time correlation functions","Duschinsky rotation","Harmonic approximation","2D resonance Raman spectroscopy","Fifth-order response functions","Finite-temperature vibronic spectra","Gaussian integrals","Duschinsky matrix"],"falsifier":"Compute the same 2DRR spectra with all sixteen Liouville pathways at the resonant carrier frequency omega_L = omega_eg and at preresonance detunings. If the twelve omitted pathways produce visible intensity or alter the antisymmetric cross-peak pattern, the paper's Duschinsky signatures are not robust under its own resonant simulation conditions.","tokens_in":18144,"feed_emoji":"🧪","tokens_out":6543,"duration_ms":65636,"temperature":0.7,"pith_summary":"The paper aims to make high-order nonlinear spectra calculable for realistic harmonic potential energy surfaces. It derives closed-form expressions for the multi-time correlation functions underlying response functions of any order, incorporating displacements, frequency changes, Duschinsky rotation between electronic states, and finite temperature without summing over vibrational eigenstates. Each correlation function reduces to a multidimensional Gaussian integral, evaluated by standard matrix operations, with polynomial scaling in the number of modes and in the response order. Applied to fifth-order 2D resonance Raman spectroscopy, the method shows that Duschinsky rotation redistributes diagonal and cross-peak intensities and creates a measurable antisymmetric spectral component, making 2DRR a sensitive probe of normal-coordinate changes.","feed_headline":"New closed forms put Duschinsky rotation into any-order spectra","feed_subtitle":"The equations include frequency shifts and finite temperature, so multidimensional Raman spectra become practical to simulate.","key_machinery":"The central object is the multi-time correlation function C(tau0,...,taun), which appears in every Liouville-space pathway of a nonlinear response function. The key reduction is to insert position-representation harmonic propagators and resolutions of identity in the ground-state normal modes, so the integrand becomes a single Gaussian in the joint (n+1)N-dimensional coordinate z. Its Hessian is the block-tridiagonal matrix D, built from the excited-state matrices a_e, b_e and the Duschinsky-rotated ground-state matrices J^T a_g J and J^T b_g J; the linear term E contains the displacement vector K. The multidimensional Gaussian integral converts C into Eq. 29/33, a ratio of determinants time","core_discovery":"The paper establishes that the nth-order vibronic correlation function C(tau0,...,taun), the building block of every Liouville-space pathway, has a closed form for arbitrary order when the electronic states are harmonic and related by the linear Duschinsky transformation q_g = J q_e + K. Writing the harmonic propagators in position representation and inserting resolutions of identity turns the trace into a Gaussian integral over (n+1)N coordinates. The result, Eq. 29 for two states and Eq. 33 for several states, is a ratio of determinants times a quadratic exponential in the displacement K, the rotated frequency matrices J^T a_g J and J^T b_g J, and the inverse of a block-tridiagonal matrix","pith_inferences":["Because the central reduction is just a Gaussian integral, the same construction should extend to non-Condon (Herzberg-Teller) dipole surfaces by inserting polynomial factors in z before integrating; this is a natural next step the paper does not carry out.","The antisymmetric cross-peak imbalance identified here could serve as a relatively background-free experimental diagnostic for normal-mode rotation, provided the four-pathway selection holds.","The polynomial scaling in response order suggests the same code could go beyond fifth order, for example to seventh-order spectroscopies, for medium-sized molecules, though the exponential growth of time-grid points with the number of delay axes remains a practical bound."],"forward_implications":["For any harmonic potential surfaces, response functions of arbitrary order can be evaluated as closed-form expressions at polynomial cost, without propagating wavepackets and without summing over vibrational eigenstates.","2DRR spectra of real molecules can be computed from standard quantum-chemistry data (frequencies, displacements, Duschinsky matrix) at finite temperature.","Duschinsky rotation changes 2DRR spectra in identifiable ways: it activates Franck-Condon-inactive modes, redistributes diagonal and cross-peak intensities, creates difference-frequency peaks, and makes the spectrum antisymmetric under exchange of the two frequency axes.","Cross peaks in 2DRR do not by themselves prove inter-mode coupling; even independent displaced harmonic oscillators produce them.","The multi-state formula allows the same machinery to be applied to processes involving more than two electronic states."],"supporting_citations":[{"why":"Supplies the nonlinear response-function and Liouville-pathway formalism, as well as the displaced-harmonic-oscillator baseline the paper generalizes.","marker":"[4]"},{"why":"Provides the practical route from quantum-chemistry data to the shift vector K and Duschinsky matrix J used to parametrize the models.","marker":"[27]"},{"why":"Defines the Duschinsky rotation of normal modes between electronic states, the physical effect the spectra are designed to probe.","marker":"[42]"},{"why":"Provides the Gaussian-wavepacket treatment that is exact for arbitrary harmonic surfaces and motivates the polynomial-scaling goal.","marker":"[54]"},{"why":"Shows how finite-temperature spectra are handled without sum-over-states, an ingredient extended here to arbitrary order.","marker":"[57]"},{"why":"Gives the squeezed-coherent-state perspective and scaling arguments that the paper's closed forms parallel.","marker":"[60]"},{"why":"Provides closed-form third-order response functions that the present arbitrary-order expressions extend.","marker":"[61]"},{"why":"Provides the naphthalene absorption benchmark used to validate the electronic-structure model.","marker":"[70]"}],"fun_headline_variants":["Closed forms for any-order spectra with Duschinsky rotation","Finite-temperature response functions go closed-form with Duschinsky","Duschinsky coupling in closed form for high-order Raman spectra","New closed formulas for finite-temp Duschinsky response functions","Any-order spectra: closed form includes Duschinsky and temp"],"cache_read_input_tokens":2688,"weakest_assumption_plain":"The 2DRR calculation keeps only four of the sixteen possible interaction histories, justified by a preresonance selection rule, yet the simulations then set the laser carrier frequency exactly on resonance; if the omitted histories contribute significantly at that resonance, the predicted Duschinsky peak asymmetries would change.","fun_headline_variants_meta":{"raw":{"variants":["Closed forms for any-order spectra with Duschinsky rotation","Finite-temperature response functions go closed-form with Duschinsky","Duschinsky coupling in closed form for high-order Raman spectra","New closed formulas for finite-temp Duschinsky response functions","Any-order spectra: closed form includes Duschinsky and temp"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000196,"raw_usage":{"total_tokens":1166,"prompt_tokens":681,"completion_tokens":485,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":425,"completion_tokens_details":{"reasoning_tokens":400}},"tokens_in":425,"tokens_out":485,"duration_ms":5211,"temperature":1.0,"reasoning_tokens":400,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-05T10:34:08.968580+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the same 2DRR spectra with all sixteen Liouville pathways at the resonant carrier frequency omega_L = omega_eg and at preresonance detunings. If the twelve omitted pathways produce visible intensity or alter the antisymmetric cross-peak pattern, the paper's Duschinsky signatures are not robust under its own resonant simulation conditions.","supporting_citations":[{"cited_title":"doi:10.1021/acs.jctc.8b00355 , issue =","cited_arxiv_id":null,"evidence_quote":"Supplies the nonlinear response-function and Liouville-pathway formalism, as well as the displaced-harmonic-oscillator baseline the paper generalizes."},{"cited_title":"Molecular Spectroscopy and Quantum Dynamics , title =","cited_arxiv_id":null,"evidence_quote":"Provides the practical route from quantum-chemistry data to the shift vector K and Duschinsky matrix J used to parametrize the models."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the Duschinsky rotation of normal modes between electronic states, the physical effect the spectra are designed to probe."},{"cited_title":"and Borrelli, Raffaele , title =","cited_arxiv_id":null,"evidence_quote":"Provides the Gaussian-wavepacket treatment that is exact for arbitrary harmonic surfaces and motivates the polynomial-scaling goal."},{"cited_title":"doi:10.1016/j.cplett.2003.09.119 , number =","cited_arxiv_id":null,"evidence_quote":"Shows how finite-temperature spectra are handled without sum-over-states, an ingredient extended here to arbitrary order."},{"cited_title":"doi:10.1021/acs.jctc.0c00079 , issue =","cited_arxiv_id":null,"evidence_quote":"Gives the squeezed-coherent-state perspective and scaling arguments that the paper's closed forms parallel."},{"cited_title":"doi:10.1063/1.3237134 , pages =","cited_arxiv_id":null,"evidence_quote":"Provides closed-form third-order response functions that the present arbitrary-order expressions extend."}],"review_version":1}