## 2 Methodology  
### Simulations  
Molecular dynamics simulations have been carried out on ternary systems composed of LiCl (where LiCl concentration ranged from 0% to 100% with a 20% interval) and (NaCl-KCl)\({}_{\text{eutectic}}\) (mole fraction ratio of NaCl:KCl is 0.506:0.494) in the temperature range from 900 to 1200 K. The simulation studies were based on the Born-Mayer-Huggins-Tosi-Fumi (BMHFT) rigid-ion (RI) interionic potential for the molten salts. The BMHFT force field reproduces the properties of molten salts well, and the formula is listed below:  
\[U_{ij}=\frac{q_{i}q_{j}}{r_{ij}}+B_{ij}exp\left(\frac{\sigma_{ij}-r_{ij}}{\rho} \right)-\frac{C_{ij}}{r_{ij}^{6}}-\frac{D_{ij}}{r_{ij}^{8}} \tag{1}\]  
where the detailed implications of these symbols are presented in nomenclature and full exposition has been provided in a previous work [13]. It should be noted that a correction factor of 0.92 was multiplied by the original calculated \(\rho\) values in the chloride melts because the parameters given by [PERSON] and [PERSON] yield shorter distances by 10-30 pm for the first peak position of the cation-anion pair correlation functions in molten ternary chlorides [6; 7]. The potential parameters used in this paper are listed in Table 1.  
The open source software, Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [14], was chosen to conduct all the simulation. A cubic box with 2000 ions (1000 cation and 1000 anion) was adopted to perform the simulation. The ion numbers of simulated compositions with a 20% interval are listed in Table 2, and were mainly discussed throughout the paper to explore the composition and temperature effect. The periodic boundary conditions method (PBC) was used to eliminate the boundary effect and keep the particle numbers constant. The cutoff radius was \(r_{c}=L/2\) where L indexed the box length. For faster calculation, the particle-particle particle-mesh solver method was employed and the precision was set as \(1\times 10^{-4}\). The [PERSON]'s equation of motion was solved by the Verlet algorithm, and the timestep during the simulation was \(\Delta\)t = 1 fs. More details of the simulated set up can be found in our previous paper [13]. For initial isothermal-isobaric ensemble (NPT) simulations, the liquid densities could be determined with the following equation:  
\[\rho_{sim}=\frac{N_{LiCI}M_{LiCI}+N_{Nucl}M_{Nucl}+N_{KCI}M_{KCI}}{v_{E}N_{A}} \tag{2}\]  
where \(N\) is the number of different alkali chlorides, \(M\) is the molar mass of each salt, \(N_{A}\) denotes the Avogadro constant, and \(v_{E}\) is the equilibrated volume of the simulation cell at the given temperature in the NPT ensemble simulations. For the structure and transport properties' collection, a random velocity according to the Gaussian distribution was assigned to the ions at the desired temperature. Then, \begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline
**Ion Pairs** & **Li-Li** & **Li-Na** & **Li-K** & **Li-Cl** & **Na-Na** & **Na-K** & **Na-Cl** & **K-K** & **K-Cl** & **Cl-Cl** \\ \hline B/eV & 0.4425 & 0.3338 & 0.3338 & 0.2905 & 0.2641 & 0.2641 & 0.2113 & 0.2641 & 0.2113 & 0.1584 \\ \(\sigma\)/Å & 1.632 & 1.986 & 2.279 & 2.401 & 2.340 & 2.633 & 2.755 & 2.926 & 3.048 & 3.170 \\ C/eV\(\cdot\)Å\({}^{6}\) & 0.04563 & 0.2189 & 0.8324 & 1.250 & 1.050 & 3.993 & 7.000 & 15.19 & 30.00 & See \\ D/eV\(\cdot\)Å\({}^{6}\) & 0.01875 & 0.09375 & 0.5313 & 1.500 & 0.5000 & 2.738 & 8.668 & 15.00 & 45.63 & below \\ \hline
**Mol Fraction of LiCl (\%)** & **0** & **20** & **40** & **60** & **80** & **100** \\ \hline \(\rho\)/Å & 0.3003 & 0.3032 & 0.3061 & 0.3091 & 0.3121 & 0.3151 \\ C\({}_{\text{ClCl-Cl}}\)/eV\(\cdot\)Å\({}^{6}\) & 75.12 & 73.98 & 72.83 & 71.69 & 70.53 & 69.38 \\ D\({}_{\text{ClCl-Cl}}\)/eV\(\cdot\)Å\({}^{6}\) & 150.9 & 148.6 & 146.3 & 144 & 141.7 & 139.4 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Potential parameters for different pairs in the ternary systems.