The role of vibronic coupling in the electronic spectroscopy of maleimide: a multi-mode and multi-state quantum dynamics study †

The first two excitation bands below 7 eV in the electronic absorption spectrum of maleimide are investigated using a model Hamiltonian including four low-lying singlet excited states within the manifold of 24 vibrational modes. The role of non-adiabatic eﬀects is studied and shines light on both the broad, inter-state coupling–dominated spectral band as well as the fine-structured, not-so-strong coupled band. Calculations have been performed using the Multiconfigurational Time-Dependent Hartree (MCTDH) wavepacket propagation method as well as its multilayer version (ML-MCTDH) using a quadratic vibronic coupling (QVC) Hamiltonian model where parameters are obtained from fitting adiabatic potential energy surfaces computed by ab initio methods. The quantum dynamics calculations provide information on the relaxation dynamics and the vibrational modes involved. Already with a low-order vibronic coupling model and only a few modes being considered, a quantitative agreement with the experimental spectrum is obtained. However, it is found that all modes need to be considered to get a full picture of the photo-excited relaxation dynamics of this molecule.


Introduction
Small heterocyclic organic compounds are ideal test systems for probing theoretical concepts in photochemistry and photophysics in order to reveal and grasp their photoactive features. This is motivated by a wide array of applications ranging from organic synthesis through drug design to polymer and biomolecular chemistry. Understanding the fundamental photochemical nature of these building blocks, with pyrrole, furan and thiophene playing a central role, facilitates the generation of new compounds with tailored properties. 1,2 In this context maleimides have gained increasing interest over the past years. In addition to the traditional use in photoinduced polymer synthesis due to their copolymerisable and self-bleaching properties, 3,4 they can be applied in a variety of research-and technology-related areas. To provide a few examples, this involves the development of new diagnostic devices 5 and antitumor agents 6 as well as the use of polyimides in the electronic industry. 7 Maleimide-based chromophores are also found in nature exhibiting an electron accepting affinity which is highly sensitive to chemical substitution. By careful substitution, the fluorescence can be tuned to a required emission frequency. 3 The impact of substitution is also seen in photoinduced reactions of several classes of maleimide derivates. 8 While maleimide is well-known for its dimerisation via intramolecular [2+2] photocycloaddition, 9 commonly observed together with photopolymerisation, N-substituted derivates avoid this reaction type and favour [5+2] photocycloadditions. 10 Derivatives containing oxygen-and thiosubstituted functional groups on the carbon double bond, known as next-generation maleimides, speed up and further expand the spectrum of accessible photocycloaddition reactions enabled by altered UV absorption mechanisms and their unique excited-state physics. 11,12 Maleimide and its N-alkyl derivates have been subject to experimental absorption spectroscopy by Seliskar et al. who identified four electronic transitions below 6.5 eV for maleimide in solution, three being allowed singlet excitations increasing in intensity. [13][14][15] These results were later confirmed by studies including a (1+1) REMPI experiment and theoretical investigation by ter Steege et al. 16 and Climent et al. 17 designating four singlet excited states of np* and pp* character: 1 1 B 1 (np*), 1 1 A 2 (np*), 1 1 B 2 (pp*) and 2 1 B 2 (pp*). Regarding the vibrational spectrum of maleimide, more theoretical [18][19][20][21] and experimental data are available including infrared, [22][23][24] Raman 23,24 and inelastic neutron scattering spectroscopy. 24 To our knowledge there are no highly-resolved vapour phase absorption spectra published so far impeding a spectral mapping with vibronic simulations.
The aim of this work is to achieve a deeper understanding of the excited-state photophysics of maleimide by studying the fine structure of the first two excitation bands in the vapour phase absorption spectrum related to the 1 1 B 2 and 2 1 B 2 excited state, respectively. Solving the time-dependent Schrödinger equation using molecular quantum dynamics simulations has been proven a solid and efficient strategy to study small aromatic and heterocyclic compounds where non-adiabatic effects are well-known to dominate the spectrum composed of vibronically-coupled electronic states. 2 By applying the vibronic coupling model of Köppel et al. 25,26 to construct an approximate Hamiltonian, the Multiconfigurational Time-Dependent Hartree (MCTDH) wavepacket propagation method [27][28][29] and its multilayer version (ML-MCTDH) [30][31][32][33][34] are used for calculating dynamical and spectral properties including a fixed number of electronic states within the manifold of vibrational modes. Although we only consider the singlet states in this work, the triplet states are also likely to play a role in the photoexcited dynamics of this molecule, and their importance will be treated in future work.

The vibronic coupling model
The vibronic coupling model has become routine in the description of vibronically coupled states by setting up an approximate Hamiltonian to cover the features required to simulate vibrationally-resolved electronic spectra. [35][36][37][38][39][40][41][42] A diabatic electronic basis is used in order to avoid singularities at points where potential energy surfaces (PESs) cross such as conical intersections. While the diabatic basis is crucial for dynamical calculations, the adiabatic picture still holds when describing the important features of the PESs following from ab initio calculations. 43,44 In the diabatic picture it is assumed that the PESs are smooth and a potential matrix W can be defined which is expressed in a Taylor series in terms of massand frequency-scaled normal mode coordinates Q. If the expansion is truncated at second order, the approach is known as the quadratic vibronic coupling (QVC) model.
The zeroth-order ground state Hamiltonian H (0) is represented by the on-diagonal nuclear kinetic energy operator T N and has the dimensionality of five, i.e. the number of electronic states being considered here. In this model, the vibrational degrees of freedom, a, enter with frequencies o a . The zerothorder diabatic potential matrix W (0) defines the main shape of the potentials. While various functions can be chosen, the simplest form is of the ground-state adiabatic potential in the harmonic approximation shifted by the relative energy E i of each excited state at the Franck-Condon (FC) point.
Up to first and second order the diabatic potential matrix elements are then defined using the expansions of eqn (3).
The linear intra-state coupling constants k (i) a are related to the gradients of the adiabatic potentials at the FC point Q 0 with respect to the nuclear coordinates, and modulate the energy gap between states along the vibrational mode corresponding to Q a . The linear inter-state coupling constants l (i,j) a provide an interaction between electronic states close in energy. The bilinear (quadratic) on-diagonal coupling constants g (i) a,b (g (i) a,a ) account for the mode-coupling Duschinsky rotation (frequency shifts). It was found that both the frequency shifts and Duschinsky rotation can become non-negligible for selected vibrational modes. The bilinear (quadratic) off-diagonal coupling constants m (i, j) a,b (m (i, j) a,a ), that provides higher order coupling between electronic states, are usually of minor importance throughout and can be reasonably neglected.
By simple group theory arguments the molecular symmetry can be used to reduce this large number of contributing parameters. Since all coupling constants can be derived as a sum of integrals with integrands belonging to a certain direct product representation, these constants will be zero if the integrands do not contain the totally symmetric irreducible representation G A of the point group of interest (here C 2v ). Thus, a set of selection rules is given depending on the electronic state symmetries For some vibrational modes, the significant anharmonicity in the potentials indicate that the harmonic approximation is a poor description and the diabatic potentials must be replaced by state-specific quartic, W q,ii , or Morse, W m,ii , potentials.
For the Morse potential, the parameter D (i) 0 defines the statespecific dissociation energy, a (i) a controls the width of the potential, and Q (i) 0,a describes the equilibrium position relative to the FC point of normal mode coordinate Q a . In the quartic potential, e (i) a is the quartic expansion coefficient. The coupling parameters for the PESs entering the model Hamiltonian are then obtained from a least-squares fitting procedure introduced in Section 3.3. It was checked that the use of Morse potentials, which are unbound, does not introduce errors into the simulations due to the wavepacket reaching the end of the grid. In all the modes for which a Morse potential was used the dissociation channel is high enough in energy not to be accessible.

The MCTDH wavepacket propagation method
Utilising the derived multi-mode and multi-state model Hamiltonian, the goal is to solve the time-dependent Schrödinger equation for the nuclear dynamics of a molecule.
The approach presented is based on the efficient MCTDH wavepacket propagation method in its mode-combined multiset framework. The wavefunction C is expanded in statedependent components C (i) indicated by the electronic state index i = 1,. . ., n e . In order to save computational effort and memory requirements, several individual modes Q km for m = 1,. . ., n k are treated together as a combined degree of freedom Q k with k = 1,. . ., p.
The state-and time-dependent multiconfigurational wavefunction C (i) can then be written in terms of a set of expansion coefficients A . . .
For numerical reasons, the SPFs are represented by timeindependent primitive basis functions w (k) l , which are a discrete variable representation (DVR) grid. For technical details on this approach, the reader is referred to the literature. [27][28][29]45,46 A limiting factor in MCTDH calculations is the length of the vector A describing the set of expansion coefficients. This length is reduced by combining several degrees of freedom in the same SPF (cf. eqn (8)). However, now high-dimensional SPFs need to be propagated and if one over combines, the propagation of the SPFs can take longer than the propagation of the A vector and the efficiency is lost. If the propagation of the multidimensional SPFs is still too long, MCTDH can be used to propagate them in a multilayer manner (ML-MCTDH). In this case, the SPFs are expanded recursively in sets of timedependent SPFs and only the last layer is represented by the discrete time-independent grid. This latter approach has been used extensively to treat large systems. [30][31][32][33][34]

Derivation of spectral properties
The simulation of the electronic excitation spectrum requires calculating the absorption cross section s, which is proportional to the squared transition matrix element of an appropriate transition operator T. The initial wavefunction C i 0 denotes the vibrational and electronic ground state which serves as reference for the vertical excitation into the excited state manifold.
In a time-dependent scheme the spectral signature follows from the Fourier transform of the autocorrelation function C. with where For a symmetric Hamiltonian and real initial wavefunction, the relation of the autocorrelation function in eqn (13) holds and the propagation length during the MCTDH calculation can be halved saving computational requirements. Fourier transform artifacts arising from a finite propagation time are reduced by introducing a cosine function, with T being the actual propagation time, to ensure that the autocorrelation function vanishes when t -T. To account for the spectral resolution, an exponential function is added to adjust the phenomenological broadening by tuning the damping time t. In the limit t -N a stick spectrum is obtained showing only discrete excitation values. Since very long propagation times are computationally inaccessible, another advantageous method has been shown to be satisfactory for many applications. 26,44 Reinserting the vibronic eigenstates into the integral form of eqn (11), the timeindependent Schrödinger equation can be recovered. When using a direct product ansatz of harmonic oscillator eigenfunctions, the unperturbed Hamiltonian T N + V now states a highly sparse algebraic eigenvalue problem that can be solved efficiently by Lanczos diagonalisation yielding the excitation energies and intensities.

Experimental absorption spectrum
The experimental absorption spectrum has been recorded in the range 3.6-6.2 eV using a PerkinElmer Lambda 365 spectrometer. The maleimide sample was heated up to 315 K in a gas cell containing fused silica windows to produce a detectable vapour pressure. The first two excitation bands occur at 4.7 eV and 5.9 eV as shown in the experimental spectra together with simulation results in Section 3.4. The first one, described in detail by Seliskar et al. in 3-methylpentane at 77 K and EPA mixed solvent at 300 K, 13 is broad and shows little vibrational fine structure fading towards higher energies. The second one has only been characterised for the N-alkyl derivates 13 indicating well-structured progressions for these type of compounds. The intensity ratio between the first and second excitation band is measured to be 1 : 15.6 which is in excellent agreement with our EOM-CCSD and CASPT2 value of 1 : 15.3 (cf. Table 2) as well as Seliskar et al.'s ratio of 1 : 13.3.

Electronic structure calculations
The diabatic PESs entering the model Hamiltonian are obtained from fitting adiabatic PESs which, in turn, result from ab initio single-point calculations along the mass-and frequency-scaled normal mode coordinates. Therefore, geometry optimisations and ground-state frequency analyses were performed at the CCSD level of theory to match experimental IR data as listed in Table 1. The mode and symmetry labels follow the energetically-ordered CCSD/cc-pVDZ results whose normal mode coordinates span the adiabatic PESs. In general, deviations get smaller when extending the basis set used. The most important normal modes are shown in Fig. 1b. It is worth mentioning that in contrast to related heterocycles such as pyrrole, 40 the N-H bond dissociation during photoexcitation, approximately described by the totally-symmetric stretching mode n 24 , is of minor importance due to a high energy barrier. Hence, the vibronic coupling model does not require further adjustments regarding this vibrational mode.
The accurate description of excited-state surfaces often requires the inclusion of static correlation in terms of multireference methods employing an active space such as CASSCF/CASPT2. 47,48 The active space CAS (12,9) has been chosen to comprise every orbital of np and pp* character. The five orbitals that mainly characterise the first five electronic states are shown in Fig. 1a. The CASSCF computations used a state average over eleven roots. Multistate CASPT2 was performed on top of the CASSCF wavefunction to account for the missing dynamical correlation. The calculations were done in absence of symmetry and no level shifts or IPEA shift were needed. Since it was not possible to produce continuous PESs for all nuclear configurations and excited states at bearable computational costs, the SA(11)-CASPT2(12,9)/cc-pVDZ and EOM-CCSD/aug-cc-pVTZ single-point results have been used to verify the applicability of the EOM-CCSD/cc-pVDZ level of theory. 47 The electronic state and symmetry labels given in Table 2 follow the energetically-ordered EOM-CCSD/cc-pVDZ results which have been chosen for the computation of the adiabatic PESs due to their appropriate accuracy-cost balance. Analysing the character of the transitions indicates that the 3s, 3p and 3d Rydberg and corresponding mixed-valence orbitals do not contribute below 7 eV, as also mentioned by Climent et al. 17 using CASPT2 calculations with varying active spaces. At all levels of the theory tested the same order of states was found, with EOM-CCSD overestimating the excitation energy by approximately 0.5 eV compared to the best CASPT2 results. It was found that only the S 3 and S 4 states are bright, with the Table 1 Experimental 22 and ab initio vibrational frequenciesñ, for the S 0 (1 1 A 1 ) ground state of maleimide in C 2v symmetry ordered by increasing energy. The CCSD/cc-pVDZ level of theory is chosen as best balance between accuracy and computational cost. Stretching vibrations are labeled by n, bending vibrations by d and out-of-plane vibrations by g followed by the most prominent atomic counterparts. The labels sym and asym define the retained and broken symmetry with respect to the s v mirror plane perpendicular to the molecular plane latter having the largest oscillator strength. These states are therefore responsible for the bands measured in the vapour phase spectrum centred at 4.72 and 5.95 eV. Cuts through the PESs at the EOM-CCSD/cc-pVDZ level of theory were made at points along all the normal mode coordinates for the chosen electronic states, denoted as S 0 -S 4 . Examples of the cuts are shown in Fig. 2. A larger state manifold comprising ten singlet excited states, S 1 -S 10 , and the first five triplet excited states, T 1 -T 5 , can be found in the ESI † employing EOM-CCSD/cc-pVDZ, along with vertical excitation energies using the SA(11)-CASPT2 (12,9), EOM-CCSD and ADC(3) methods. Coupling between the S 4 and S 5 excited state is vanishingly small throughout owing to the big energy gap at the FC region. Thus, ignoring the triplet states, a model Hamiltonian selecting the first five singlet electronic states is expected to cover most of the vibronic features necessary for the description of the absorption spectrum below 7 eV.
The characterisation of excited-state PESs along selected normal mode coordinates can be divided into two classes of vibrational modes that contribute to the model Hamiltonian (cf. Fig. 2). Mode n 18 is representative for symmetric stretching vibrations of a 1 symmetry that provide non-zero intra-state coupling and show surface crossings. Mode n 1 class vibrations on the other hand correspond to either symmetric out-of-plane or asymmetric bending vibrations of b 1 or b 2 symmetry. These potentials are symmetric with respect to the FC point and these modes can provide inter-state and higher-order coupling (cf. Section 3.3). The double well of the S 2 surface along Q 1 in Fig. 2 is due to strong inter-state non-adiabatic coupling with the S 3 state. An overview of the PESs of all 24 normal modes can be found in the ESI. † The CCSD geometry optimisation and frequency analysis, as well as excited state calculations employing the EOM-CCSD method, have been performed using Gaussian16 49 and Q-Chem5.3. 50 The CASSCF/CASPT2 calculations were performed with the program Molcas8.2. 51

Fitting the potential energy surfaces
The coupling parameters which need to be determined are included in the model potentials W entering the electronic Hamiltonian in the diabatic basis. Since diagonalisation of the corresponding matrix provides the adiabatic energies, the coupling parameters can be determined by minimising the root-mean-square deviation of these energies with respect to the computed single-point energies at different nuclear configurations. This was done using the VCHam programs, which are a part of the quantum dynamics software Quantics. 52,53 The actual quantity to be optimised (eqn (16)) is the functional D with s being here the number of electronic states and n s the number of points along each adiabatic PES considered. A weight function O s i is added where V s 0 corresponds to the sth adiabatic potential at the FC point.   13 , ab initio vertical singlet excitation energies, E, and oscillator strengths, f, in parentheses for maleimide in C 2v symmetry. The EOM-CCSD/cc-pVDZ level of theory is chosen for calculating points on the adiabatic PESs, while the SA(11)-CASPT2(12,9)/cc-pVDZ and EOM-CCSD/aug-cc-pVTZ levels of theory are used for verifying obtained results. The first five electronic states S 0 -S 4 are considered in the model Hamiltonian. The state characters follow from the largest CI vectors of the corresponding CASPT2 calculation State Symmetry EOM-CCSD/cc-pVDZ EOM-CCSD/aug-cc-pVTZ SA(11)-CASPT2(12,9)/cc-pVDZ Experiment The optimisation is performed order-by-order in the Taylor expansion with on-diagonal elements being calculated before off-diagonal contributions. It is known that linear intra-and inter-state coupling constants k (i) a and l (i,j) a , listed in Tables 3  and 4, are responsible for the main features of the spectra while the quadratic and bilinear on-diagonal coupling constants g (i) a,a and g (i) a,b only have a small impact for most vibrational modes. The most important normal modes based on first-order parameters are shown in Fig. 1b. The full Hamiltonian is provided in the ESI. † Since the intra-state coupling constants k (i) a correspond to the gradients of the adiabatic potentials at the FC point, large non-zero values weighted by the corresponding frequencies, o a , indicate modes that will significantly influence the dynamics following the vertical excitation to the S 3 (1 1 B 2 ) and S 4 (2 1 B 2 ) singlet excited states. Only totally-symmetric a 1 modes exhibit non-zero gradients at the FC point, as illustrated for mode n 18 in Fig. 2 (cf. eqn (4)). Concerning the vertical excitation to the S 3 (1 1 B 2 ) state, modes n 7 , n 18 and n 21 are of special importance as they exhibit surface crossings with the S 2 (1 1 A 2 ) state which provide a barrierless pathway for the relaxation. Although they do not lead directly to any crossings, n 7 , n 18 , n 19 and n 21 are the most important vibrations for excitation to the S 4 (2 1 B 2 ) state.
Non-zero inter-state coupling constants l (i,j) a represent the first-order off-diagonal elements of the diabatic potential matrix meaning that these are the most important nonadiabatic coupling contributions with respect to the shorttime dynamics. The larger the frequency-weighted values are for electronic states close in energy, the bigger is the impact on the vibronic spectrum. Significant coupling is obtained between the S 3 (1 1 B 2 ) and lower lying S 2 (1 1 A 2 ) and S 1 (1 1 B 1 ) excited states. Specifically the vibrational modes n 1 (b 1 ), n 2 (a 2 ), n 4 (b 1 ), n 9 (a 2 ) and n 20 (b 2 ) dominate the coupling between these energetically close states. The excited state S 4 (2 1 B 2 ) is only coupled weakly to the low-lying states due to the increased energy gap around the FC region. As a consequence, first-order non-adiabatic effects are small in the second excitation band leading to a clear fine structure in the spectrum. Nevertheless, significant couplings are present with the S 3 (1 1 B 2 ) state along Fig. 2 Adiabatic single-point energies, E, (dots) of maleimide calculated at the EOM-CCSD/cc-pVDZ level of theory for the first five electronic states shown relative to the ground state energy at the FC point, E 0 (S 0 ), as a function of mass-and frequency-scaled normal mode coordinates Q of mode n 1 (left) and n 18 (right). The adiabatic potential energy curves (solid lines) are obtained from the fitted vibronic coupling model for that particular normal mode. Both the single-point energies as well as the fitted curves have been shifted to the SA(11)-CASPT2(12,9)/cc-pVDZ excitation energies. The ground state S 0 is shown in black, S 1 in yellow, S 2 in green, S 3 in red and S 4 in blue with symmetry labels referring to maleimide in C 2v symmetry at the FC point. the n 3 and n 7 vibrations. Note that these are totally symmetric vibrations as the S 3 and S 4 states have the same symmetry and these modes have significant k parameters as well. Second-order on-diagonal contributions of quadratic and bilinear character have also been included. The quadratic constants g (i) a,a affect the spectrum by a state-specific change of frequency for mode a as shown in eqn (17).
The bilinear constants g (i) a,b couple the modes a and b and describe Duschinsky rotation. This effect is the mixing of the vibrational modes of one electronic state in the normal mode basis of the ground electronic state. These were all found to be small. By inspection, modes n 1 , n 4 and n 22 are not harmonic and were described instead by a quartic oscillator diabatic potential, as shown in eqn (5). Modes n 11 , n 19 , n 21 , n 23 and n 24 were described by the Morse potential of eqn (6). The parameters obtained for these potentials from the least-squares fit are given in the ESI. † It was found that minor adjustments had to be made to some parameters to achieve a good agreement with the experimental spectrum, in particular to get the fine structure of the second spectral band. These adjusted parameters are presented in Table 5. In addition, the CASPT2 energies were used to adjust the position of the EOM-CCSD potential energy surfaces at the FC point.

Quantum dynamics calculations
Within the quantum dynamics simulations employing the MCTDH method, a wavepacket is vertically excited to the singlet states S 3 (1 1 B 1 ) and S 4 (2 1 B 1 ) from which the absorption spectra are obtained by evaluating the autocorrelation function, i.e. a constant transition dipole moment is assumed. The spectra are then normalised to have the same maximum intensity. Two models have been set up for the MCTDH simulations containing six and twelve modes, as well as a full 24-mode model using the ML-MCTDH approach. The 6-mode model for the simulations starting on the S 4 state includes all of the totally symmetric modes, which provide the intra-state coupling and coupling to the S 3 state. For the dynamics starting on the S 3 state, the n 3 and n 21 modes that have the weakest inter-state coupling are replaced by n 1 and n 4 to include the coupling to the S 2 state.
The 12-mode model Hamiltonian is the same for the simulations starting on both states and includes all of the modes with strong coupling terms listed in Tables 3 and 4, except n 10 that should have the smallest effect on the dynamics as it only weakly couples states S 2 and S 4 . Details concerning the mode combination, the number of SPFs and the DVR grid for converged spectra and state populations up to 200 fs can be found in Table 6 and Fig. 3. All calculations were done using the Quantics Package. 52,53 3.4.1 The first 1 1 B 2 ' 1 1 A 1 excitation band. Fig. 4 presents the experimental vapour phase absorption spectrum for the first band centred at 4.72 eV, along with the simulated 6-mode, 12-mode and 24-mode model results starting in the S 3 state. The experimental spectrum shows the same vibrational progressions recorded by Seliskar et al. 13 although being more pronounced and shifted due to the absence of a solvent. The 6-mode model comprises the most important intra-state constants of a 1 modes n 7 , n 11 , n 18 and n 19 , along with the strongest coupling to the S 2 state through the b 1 modes n 1 and n 4 . However, the broad nature of this excitation band is not completely described by this model, as shown by the spectra produced with a damping time of  Table 5 Selected parameters which require an adjustment compared to those obtained from the fitting procedure with the objective to reproduce the experimental spectra best. The first three values refer to the anharmonic Morse potentials and the last three values to the coupling parameters arising from the Taylor expansion in the QVC model. All potentials have been shifted to the SA(11)-CASPT2 (12,9)

/eV
À0.061 À0.100 g (4) 18,18 /eV À0.014 À0.020 either 30 or 200 fs: the former is required to wash out the structure, indicating missing couplings. The resulting spectra with 30 fs damping is in good agreement with the experimental one, with the same width and similar weak structure. With a damping time of 200 fs the vibronic lines are enhanced. It is impossible to resolve them further. The spectrum is thus found to be a complicated set of multi-mode lines and it is difficult to assign the structure to particular vibrations as the bands overlap. The main shape and width, however, is due to the strong excitation of the n 7 , n 11 and n 18 modes. These are responsible for the strong progression of peaks -the first four of which are at approximately 0.08, 0.12, 0.17, 0.20 and 0.25 eV above the origin (marked with an asterisk) and thus due to 7 1 0 , 11 1 0 , 7 1 0 + 18 1 0 , 7 1 0 11 1 0 and 7 1 0 18 1 0 vibronic excitations. There is also clear structure below the S 3 origin due to excitations involving the S 2 state.
The 12-mode model adds in the remaining a 1 modes n 3 and n 21 , and includes the inter-state coupling to the S 1 state, either directly to S 3 or indirectly via S 2 , by adding the modes n 2 , and n 20 , and the coupling to S 0 by modes n 9 and n 13 . A damping time of 200 fs now produces the broad band with little structure. The full 24-mode model then masks the structure further and reproduces the experimental spectra very well. Going from the 12-mode to 24-mode models does not increase the interstate coupling appreciably, and so this extra damping of the spectrum must be due to intra-state dephasing in the full set of vibrations.
The diabatic population transfer for the three models between the excited states S 1 -S 4 is shown in Fig. 5. While adiabatic populations refer to electronic state surfaces provided by the clamped-nucleus Hamiltonian, the diabatic populations are based on PESs for a certain electronic configuration. 44 For example the S 3 (1 1 B 2 ) diabatic state is the p 2 p* configuration, and the S 2 (1 1 A 2 ) the n 1 p*, etc. These populations can thus be related to the electronic structure of the state at the FC point and are obtained directly from the wavepacket dynamics simulations due to the construction of the model.
Considering the 6-mode model first, an ultrafast decay of the initially populated S 3 (1 1 B 2 ) state is obtained due to the coupling with the adjacent S 2 (1 1 A 2 ) state. This initial transfer of population can be associated with the conical intersections between these states close to the FC point. At around 100 fs the populations of those states becomes approximately saturated with a transfer of 60% to the S 2 state, while the S 4 (2 1 B 2 ) makes a negligibly small contribution. The S 1 and S 0 states are not coupled and do not receive any population. The populations in states S 2 and S 3 are seen to coherently oscillate with two time period of approximately 25 and 50 fs, which correlates nicely with the time periods of n 7 and n 18 which are approximately 25 and 50 fs, respectively.
The calculations performed with more modes, i.e. the 12and 24-mode models, show that adding the coupling to the lower lying states dramatically changes the population dynamics and now the population is transferred mostly to the S 1 state. With 12-modes, after 100 fs there is approximately a ratio of  [2,9,4,14,9] HO (16,21,11,31) (n 2 ,n 19 ,n 20 ,n 21 ) [2,9,3,13,8] HO (11,16,11,16) Fig . 3 Layer structure of the ML-MCTDH simulations including all 24 normal modes. The number in the circles (grey) denotes the node number in the layered structure, while within nodes the number n of SPFs is noted next to the link lines. From the second layer onwards two values are given representing the simulation for the first and second excitation band, respectively. In the last layer, comprising the normal modes of one ''particle'', the number N of primitive functions (blue) is used to represent the grid in each normal mode. The maximum depth of the tree is five layers and the first one separates the 24 vibrational coordinates from the discrete electronic degree of freedom (red) in the single-set formulation. 31 66 : 12 : 20 : 2 in the S 1 , S 2 , S 3 and S 4 states. A negligible transfer to S 0 is found. While the initial transfer goes from S 3 to S 2 , a coherent transfer between S 3 and S 1 with a time period of approximately 50 fs also takes place. With all 24 modes included, a very similar population dynamics is seen, but the final ratios are now approximately 72 : 12 : 12 : 3 for S 1 : S 2 : S 3 : S 4 , and less than 1% ends up in the ground-state. Thus, the many small couplings play a crucial role in the population transfer. This can be explained by the inclusion of the vibrational modes n 2 (a 2 ), n 9 (a 2 ) and n 20 (b 2 ) which mediate the population transfer into S 1 , either directly from S 3 or via S 2 . Hence, even though the 6-mode model can qualitatively reproduce the experimental spectrum, a higher-dimensional Hamiltonian is essential for describing the full dynamical behavior.
3.4.2 The second 2 1 B 2 ' 1 1 A 1 excitation band. Fig. 6 presents the experimental absorption spectrum of the second intense band in maleimide along with the simulated 6-mode, 12-mode and full 24-mode model results. In contrast to the first excitation band, the spectrum corresponding to the energetically more isolated S 4 (2 1 B 2 ) state is sparser and has more structure. Hence, the main contributions should arise from totally-symmetric a 1 modes with strong intra-state coupling. Hot bands at the low energy end of the experimental spectrum have been identified by tracking their intensity ratios when heating the sample to different temperatures. The simulated spectra have all been broadened with a damping time of 200 fs.
As the S 3 and S 4 states have the same symmetry, the totally symmetric modes in the 6-mode model include the linear coupling between these states. As a result, the 6-mode model comprising the a 1 modes n 3 , n 7 , n 11 , n 18 , n 19 and n 21 is capable of producing results of quantitative accuracy for the energy regime below 6.0 eV. The low-energy regime is mostly shaped by fundamental vibronic excitations as labeled in Fig. 6. A more detailed peak assignment is given in the ESI. † The highenergy regime of the spectrum, however, shows a more complex  13 The simulated spectrum results are shown for the vertical excitation to the S 3 (1 1 B 2 ) diabatic state using a 6-mode (dotted line), 12-mode (dashed line) and 24-mode model (solid line). The simulated spectra have been shifted to overlap the experimental band by approximately the zero-point energy of the ground-state. A damping time of t = 30 fs (red) has been used for the 6-mode model to account for the phenomenological broadening. The * indicates the origin of the vibrational progressions found in the 6-mode model. The vibronic fine structure for the 6-mode model is enhanced using a damping time of t = 200 fs (light red) and the stick spectrum in the limit t -N (grey) is obtained from direct Lanczos diagonalisation in 1000 steps. The 12-mode and 24-mode spectra are broadened using only a 200 fs damping time. The spectra are shifted vertically for reasons of visibility. state (red) into the S 1 (1 1 A 2 ) state (yellow), the S 2 (1 1 B 1 ) state (green) and the S 4 (2 1 B 2 ) state (blue). The S 0 (1 1 A 1 ) ground state is omitted due to negligible diabatic population. The initial 50 fs region (grey area) is shown enlarged below. structure of broader nature than the experimental one, which can be explained by the errors in the Hamiltonian parameters and would be improved by further adjustments.
Adding the extra modes in the 12-mode and 24-mode models does not affect the main features of the spectrum, except for increasing the intensity in the middle of the band. It should be mentioned that heating or cooling the sample would introduce a temperature-dependent broadening of the experimental peak widths as well as producing a tail towards higher energies additionally affecting the absolute intensities.
The diabatic population transfer following the excitation to the S 4 (2 1 B 2 ) state is depicted in Fig. 7. It is clear to see that by just considering the dominant intra-state coupling, the wavepacket preferably remains on the initial excited state and only about 5% transfer to the S 3 occurs. Adding in the main couplings to the lower states in the 12-mode model does not appreciably change this behaviour, except that now some population is transferred to S 1 , with a ratio of approximately 2 : 0 : 9 : 89 after 200 fs for S 1 : S 2 : S 3 : S 4 .
There is, however, a striking change in the population dynamics when all the modes are included. There is now a steady transfer of population to the lower states, and after 200 fs the ratio of populations is approximately 25 : 10 : 20 : 45 in the S 1 , S 2 , S 3 and S 4 respectively. This is a huge change due to the addition of many small couplings between the states (listed in ESI †). Looking in more detail at the initial population dynamics, it can be seen that at the start the additional modes provide a slightly larger transfer to S 3 and S 2 states. After 10 fs, this leads to a steady flow of population into S 1 and S 2 , while the S 3 population starts to plateau.

Conclusions
This work constructs an approximate model Hamiltonian in order to study quantum dynamical aspects of the first two excitation bands of maleimide's electronic absorption spectrum. The experimental spectrum is shown for our vapour phase result at 315 K (black). The simulated spectrum results from the vertical excitation to the S 4 (2 1 B 2 ) diabatic state using a 6-mode (dotted line), 12-mode (dashed line) and 24-mode model (solid line). The simulated spectra have been shifted to align with the experimental band by subtracting approximately the ground-state zero-point energy. A damping time of t = 200 fs (blue) has been chosen to account for the phenomenological broadening. The stick spectrum in the limit t -N (grey) is obtained from direct Lanczos diagonalisation in 1000 steps for the 6-mode model. The simulated spectra are shifted for reasons of visibility. The peak assignment of the most important vibronic excitations in the experimental spectrum follows from the 6-mode model (cf. ESI † for details). ) state (blue) into the S 1 (1 1 A 2 ) state (yellow), the S 2 (1 1 B 1 ) state (green) and the S 3 (1 1 B 2 ) state (red). The S 0 (1 1 A 1 ) ground state is omitted due to negligible diabatic population. The 50 fs region (grey area) is shown enlarged below.
In short, the quadratic vibronic coupling model is applied together with the MCTDH wavepacket propagation method. Adiabatic potential energy surfaces, from which the diabatic electronic basis is obtained, were calculated at the EOM-CCSD level of theory and verified with CASPT2. Only a small number of parameters required adjustment to provide good agreement between the model and experiment.
New experimental spectra have also been recorded. The good agreement between the experimental absorption spectra and the simulated spectra from quantum dynamics simulations following vibronic excitation can be regarded as a validation of the simple model Hamiltonian. The two excitation bands studied have been found to be of fundamental different character. The 4.72 eV transition is dominated by inter-state coupling, and is broad and structureless. In contrast, the 5.95 eV transition has detailed structure owing to weaker inter-state couplings.
The second band is clearly identifiable as being due to excitation to the S 4 (2 1 B 2 ) state, and is well reproduced by a 6-mode model that includes all the totally symmetric vibrations that provide both the intra-and inter-state couplings to the lower S 3 (1 1 B 2 ) state. Little change is made to the spectrum by adding more modes and it was possible to assign the spectral peaks. The first band is identified as being due to excitation to the S 3 (1 1 B 2 ) state. Both bands are pp* transitions. In contrast to the S 4 band, the full 24-mode model was required to simulate the broad, fairly featureless spectrum. Lower dimensional calculations could, however, give insight into the vibrational modes being excited and the origin of the structure in the spectrum was found to be multi-mode in nature, with two dominant vibrations.
The population dynamics of the two bands also shows a different behaviour. Relaxation from the S 3 state is fast and goes directly to both the S 2 and S 1 states. Negligible population returns to the ground state. Both the 12-mode and 24-mode models show the same dynamics, with over 70% of the population transferred to the S 1 state after 100 fs. In contrast, relaxation from the S 4 state in both the 6-mode and 12-mode models is seen to be very small. Only once all of the modes with many small inter-state couplings are added, an appreciable relaxation is seen (55% after 200 fs). Most of the population ends up in the S 3 and S 1 states.
In conclusion, the four lowest singlet excited states of maleimide provide a strongly coupled manifold. Even from the energetically more separated S 4 state, which leads to a highly structured spectral band, ultrafast vibronic relaxation is seen providing all the vibrational modes are taken into account. Interestingly, after vertical excitation to the S 3 state there is a highly efficient pathway directly to the S 1 state, but after vertical excitation to the S 4 the state population is left trapped in the S 3 as well as relaxing down to the S 1 state.

Conflicts of interest
There are no conflicts to declare.