As pointed out in Sec. 2.3, the dipoles, including the permanent dipoles and transition dipoles, and their derivatives are necessary ingredients to perform polariton dynamics simulations. However, for realistic molecular systems, the derivatives of dipoles are not readily available for most of the commonly-used excitedstate electronic structure methods, such as CASSCF and TD-DFT. Therefore, we circumvent this technical difficulty by employing the machine learning techniques in this work to get the analytical expressions of dipoles in terms of the molecular geometry µ(R) = f (R). After that, we have full access to the analytical expression of derivatives of dipoles using ∇µ(R) = ∇f (R). Below, we will briefly discuss the machine learning strategy we used to parametrize dipoles, and the details can be found in the previous work on performing on-the-fly non-adiabatic dynamics with the machine-learned potential energy surfaces (PESs). To obtain the relation between the dipoles and the molecular geometry, we first need to define a proper molecular descriptor to represent the molecular geometry. One simple molecular descriptor is the Coulomb matrix M, with the matrix elements defined as follows (48) where Z k is the atomic number of atom k, (x k , y k , z k ) is the Cartesian coordinates of atom k. Note that the diagonal elements of the Coulomb matrix are not dependent on molecular configuration, and thus it is safe to remove them in the construction of molecular descriptor. In addition, the Coulomb matrix is symmetric, thus we only need to consider the off-diagonal elements of the upper triangle part. All of these used off-diagonal elements define a vector m that contains N a × (N a -1)/2 elements (where N a is the number of atoms of a single molecule) and are used as input in the training and prediction process. The Coulomb matrix provides a simple and effective representation of molecular geometry, which takes both element types and internal distances into account and has shown great advantages in building machine learning models as demonstrated in the previous work. In this work, we employ the KRR method in the machine learning process. KRR is one of the most popular supervised learning approaches and has been used for the prediction of molecular properties in several studies. In the KRR approach, the molecular property f (m), such as the molecular dipole considered in this work, is estimated by a function of the nuclear configurations represented by the molecular descriptor m. The model assumes that the property can be calculated as the follows