A method is presented for the unbiased numerical computation of two-particle response functions of correlated electron materials via a solution of the dynamical mean-field equations in the presence of a perturbing field. The power of the method is demonstrated via a computation of the Raman B1g and B2g scattering intensities of the two dimensional Hubbard model, in parameter regimes believed to be relevant to high-temperature superconductivity. The theory reproduces the 'twomagnon' peak characteristic of the Raman intensity of the insulating parent compounds of high-Tc copper oxide superconductors and shows how it evolves to a quasiparticle response as carriers are added. The method can be applied in any situation where a solution of the equilibrium dynamical mean-field equations is feasible. [3,4], has transformed our understanding of correlated electron physics. In the particular case of the two-dimensional Hubbard model, believed [5] to represent the physics of high temperature superconductors, the method has provided new insights into the correlation-driven ('Mott') metal-insulator transition [6,7], the pseudogap regime that separates the Mott insulator from the Fermi liquid metal [7][8][9][10] In this paper we present a new method for determining the two-particle response in cluster dynamical mean-field theory and demonstrate its effectiveness via a computation of the doping dependence of the Raman scattering amplitude of the two-dimensional Hubbard model from the Mott insulating to the Fermi liquid regime. Raman spectroscopy has been of fundamental importance to high temperature superconductivity [24] but the theoretical description in terms of an underlying Hubbard model involves vertex corrections in an essential way and has not been systematically studied.To introduce our method we recall salient features of the theory of linear response [25], defined quantum mechanically as the leading-order difference of the expectation value of an operatorR in the presence and absence of a probe field P . This is given by R P = Tr R G P − Tr R G eq = χ RP P + O(P 2 ) with(1) Here G eq = (G −1 0 − Σ eq ) −1 is the equilibrium (P = 0) Green function, related in the standard way to a bare Green function G 0 and a self energy Σ. We have omitted a possible term arising from explicit dependence of R on P ; this gives rise to the 'diamagnetic' term in the optical conductivity but is not otherwise relevant. The first term in Eq. 1 gives the 'bubble term', computable from knowledge of the one-electron Green function and the bare vertex δG −1 0 /δP ; the second term, arising from changes in the many-body physics due to the perturbation, is the vertex correction of interest here.For wide classes of strongly correlated materials, neither perturbative diagrammatic expansions about a mean-field solution nor partial (e.g RPA or GW) resummations suffice; a fully nonperturbative treatment is required. For the one electron Green function, cluster dynamical mean-field theory [2] provides such a treatment. In this theory the electron ...