Our analytic scheme uses as a starting point that linear response functions described by perturbations a and b can be formulated as perturbation strength ε i (i = a, b,…) derivatives of a quasienergy Lagrangian gradient, expressed in a density-matrix (D) formulation as48 1where the derivative is evaluated at zero perturbation strength and where higher-order response functions can be found by further differentiation of eqn (1). A tilde is used to represent a quantity considered at an arbitrary perturbation strength, and the absence of a tilde denotes evaluation at zero perturbation strength. The quasienergy Lagrangian L a is given by2where we have introduced the atomic orbital (AO) overlap matrix S as3Sμν = χ̃μ|χ̃ν,where χ̃ is an atomic orbital, and where the energy- and frequency-weighted Fock matrix W is defined as4where the generalized Kohn–Sham Fock matrix is given by5 We also introduced the generalized energy as6 7 In eqn (5)–(7), we introduced the half-time-differentiated overlap matrix T, the one-electron matrix h, the external field operator t and the two-electron matrix G γ with γ-fractional exchange as8 9 10 11and also the exchange–correlation contributions F xc and xc[ρ̃(D)] in addition to a nuclear potential operator h nuc. Here and throughout the paper, atomic units are used unless otherwise stated. Molecular properties characterized by a perturbation tuple abc… can therefore be formulated as derivatives of the quasienergy Lagrangian gradient as12 13 14where we have introduced a short-hand notation for differentiation and tracing by15and16respectively. This theory is sufficient to define any response function using the so-called n + 1 rule formulation,55 where the calculation of a response property of order n + 1 requires the calculation of the density matrix perturbed to order n. However, other formulations placing other conditions on which perturbed density (and Fock) matrices must be calculated are possible.55 Let us represent the idempotency of the density matrix and the time-dependent self-consistent field (TDSCF) conditions as the matrices Y and Z, respectively, so that17Y = DSD – D,and18where the notation19[M]⊖ = M – M†,and20[M]⊕ = M + M†,has been introduced, and where adjungation is defined to happen before time differentiation. It can be shown that the ansatz21λ̃a = [DaSD]⊖,for the multiplier λ̃ a for Y leads to the definition of the multiplier ζ̃ a for Z as22 It is then possible to make a general expression for the quasi-energy Lagrangian for the calculation of response properties as23where the values of k and n in the various terms denote, with minor variations, to which orders perturbed Fock and density matrices must be calculated in order to evaluate this expression: the value of k determines to which order the perturbed matrices must be calculated for perturbation tuples involving perturbation a, whereas the value of n determines the same for perturbation tuples not involving perturbation a. We have that k + n = N – 1, where N is the order of the property considered, and k must be chosen as an integer in the interval k ∈ [0,(N – 1)/2], where (N – 1)/2 is rounded down for even N. In this work, we do not discuss how the necessary perturbed Fock and density matrices can be calculated, as it is described in detail in ref. 48. We remark, however, that since the calculation of high-order properties requires solving linear response equation systems and since this part of the calculation is computationally expensive, a judicious (k,n) rule choice may give a significant reduction in the number of such systems to be solved, both compared to other rule choices and to numerical differentiation schemes. For instance, for the calculation of cubic force constants, the choice (k,n) = (1,1) makes it necessary to solve M systems, where M is the number of geometrical coordinates, whereas (k,n) = (0,2) results in M 2 such systems. Similarly, a scheme where an analytically calculated molecular Hessian is differentiated numerically by nuclear displacements results in the number of such systems being of the order of M 2. Similar savings can be achieved for other properties.