Modeling Coloration of a Radiochromic Film with Molecular Dynamics-Coupled Finite Element Method

: Radiochromic films change color upon exposures to radiation doses as a result of solid-state polymerization (SSP). Commercially available radiochromic films are primarily designed for, and have become widely used in, clinical X-ray dosimetry. However, many intriguing properties of radiochromic films are not yet fully understood. The present work aimed at developing a theoretical model at both atomic and macroscopic scales to provide a platform for future works to understand these intriguing properties. Despite the fact that radiochromic ﬁlms were primarily designed for clinical X-ray dosimetry, dose-response curves for the Gafchromic EBT3 ﬁlm obtained for ultraviolet (UV) radiation were employed to develop our model in order to avoid complications of ionization, non-uniform energy deposition, as well as dispersed doses caused by secondary electrons set in motion by the indirectly ionizing X-ray photons, which might introduce added uncertainties to the model and overshadow the basic SSP processes. The active layer in the EBT3 ﬁlm consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using molecular dynamics (MD). The degrees of SSP in the atomic scale upon different UV exposures were obtained to determine the absorption coefﬁcients of the active layer, which were then input into the ﬁnite element method (FEM). The classical steady-state Helmholtz equation was engaged to model the reﬂection from the active layer using the FEM technique. The multifrontal massively parallel sparse direct solver (MUMPS) was employed to solve the present numerical problem. Very good agreement between experimentally and theoretically obtained coloration in terms of net reﬂective optical density was achieved for different UV exposures. In particular, for UV exposures larger than ~40 J/cm 2 , the reﬂected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization, and thus the lower efﬁciency for further bond formation between DA monomer strands. assumed to did conclusions, as explained below). The density ρ the polymer MD computations was used to determine the absorption coefﬁcient α for i.e., µ ( µ α ) m . The variation in the absorption coefﬁcient versus different UV exposures reﬂected the atomic-scale dynamics


Introduction
Radiochromic films change colors upon exposures to radiation doses as a result of solid-state polymerization (SSP) [1], where energy is transferred to the receptive part of the monomers [2]. Commercially available radiochromic films are primarily designed for, and have become widely used in, clinical X-ray dosimetry [2]. Application of radiochromic films for clinical X-ray dosimetry has, in essence, relied on dose-response curves (response in terms of film coloration) obtained through calibration. However, many intriguing properties of radiochromic films are not yet fully understood, including their responses to fluorescent light [3,4] and ultraviolet (UV) radiation [5][6][7][8]; the energy dependence of their responses [9][10][11][12][13]; the effects of polarization and scanning orientation [14][15][16]; and the effects of mechanical pressure [17] and auto-development [18], etc. In order to understand these intriguing properties, development of a theoretical model at both atomic and macroscopic scales is pertinent. To the best of our knowledge, there are at present no theoretical models in the literature that can reproduce experimentally obtained dose-response curves of radiochromic films. The present work was devoted to developing such a model. Despite the fact that radiochromic films were primarily designed for clinical X-ray dosimetry and that most dose-response curves were understandably obtained for X-rays, development of a theoretical model in the atomic scale using X-ray data would involve added complications of ionization, non-uniform energy deposition, as well as dispersed doses [19] caused by secondary electrons set in motion by the indirectly ionizing X-ray photons. Such complications will inevitably introduce added uncertainties to the model and might overshadow the basic SSP processes. Fortunately, dose-response curves for the Gafchromic EBT3 film have also been obtained for UV radiation [8], for which complications caused by secondary electrons described above are insignificant. These data would provide reference information for the development of our model. A future paper will cover further investigations that add processes involving secondary electrons created by X-ray photons into the model developed here.
The present task was divided into two parts, namely, (1) the atomic scale model and (2) the macroscopic scale model, in which molecular dynamics (MD) and finite element method (FEM) simulation techniques were coupled. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers [20] (hereafter referred to as DA monomers), which were modelled using the MD technique. The degrees of SSP in the atomic scale upon different UV exposures were obtained in order to determine the absorption coefficients of the active layer, which were then input into the FEM model. The classical steady-state Helmholtz equation was engaged to model the re-emittance (i.e., reflection) from and transmittance through the active layer using the FEM technique. The multifrontal massively parallel sparse direct solver (MUMPS) [21] was employed to solve the present numerical problem.

Atomic Scale Model
To understand the degree of polymerization in an EBT3 film, we employed MD computations to model the behavior of DA monomers using the AIREBO (Adaptive Intermolecular Reactive Empirical Bond Order) potential [22]. The DA monomers were modeled in a rectangular simulation box with lengths of 800 Å in the x direction, 230 Å in the y direction, and 150 Å in the z direction, and the carboxylic chain (-COOH) was neglected. The number of oxygen atoms in the DA monomer was insignificant compared to those in the long hydrocarbon chain. The dominant bonding sites would be located along the simulated chain and the contribution of the -COOH group would be negligible. The attraction forces among the DA monomers were governed by attraction forces among atoms in the chains, which were inversely proportional to separations among the atoms. The simulation box consisted of~10 6 (carbon + hydrogen) atoms in total, which accounted for~15,000 DA monomer strands. A section of the simulation box is shown schematically in Figure 1.
In the AIREBO potential used in the present work, each pair of covalently-bonded atoms interacted via a potential given by where V ij R and V ij A were pairwise potentials for repulsion (R) and for attraction (A), respectively; and where the monomer chain was considered to consist of carbon (i) and hydrogen (j) atoms. The AIREBO potential consisted of three terms: where E ij LJ represented the intermolecular interaction potential, which accounted for intermolecular repulsion and dispersion defined using the Lennard-Jones (LJ) 12-6 potential: where was the potential well depth; σ was the distance at which the intermolecular potential between the two particles (i) and (j) was essentially zero; and r represented the distance between two attracting particles. Due to the presence of the torsional interaction potential in the AIREBO potential for dihedral angles, this potential would be highly suitable for long-chain hydrocarbon system such as DA monomers. The torsional potential was expressed as with ω as the dihedral angle. The energy minimization of the system was performed for 10 3 fs for DA monomers using the constant-energy, constant-volume ensemble (NVE) at 295 K. The structure was then relaxed using Langevin thermostat at 295 K for 10 3 fs, which enabled the atoms and the overall monomer chains to relax. Furthermore, the relaxed DA strands were exposed to specific UV exposures; this exposure was supplied in the form of energy to the atoms in every chain in small increments over a relatively long relaxation time, mainly to avoid excessive energy supply to the atoms that might break them from the main DA chain. After the energy supply step to the chains was completed, the DA monomers were relaxed under constant-pressure, constant-temperature ensemble (NPT) for~2000 fs, and the final density of the polymerized system was obtained. The mass energy-absorption coefficient ((µ α ) m ) for the present polymeric system was assumed to be 100 or 10 cm 2 /g (which actually did not affect the final results and conclusions, as explained below). The density ρ of the polymer system obtained from MD computations was used to determine the absorption coefficient µ α for different UV exposures, i.e., µ α = ρ × (µ α ) m . The variation in the absorption coefficient versus different UV exposures reflected the atomic-scale dynamics of the DA monomers.
dose-response curves were understandably obtained for X-rays, development of a theoretical model in the atomic scale using X-ray data would involve added complications of ionization, non-uniform energy deposition, as well as dispersed doses [19] caused by secondary electrons set in motion by the indirectly ionizing X-ray photons. Such complications will inevitably introduce added uncertainties to the model and might overshadow the basic SSP processes. Fortunately, dose-response curves for the Gafchromic EBT3 film have also been obtained for UV radiation [8], for which complications caused by secondary electrons described above are insignificant. These data would provide reference information for the development of our model. A future paper will cover further investigations that add processes involving secondary electrons created by X-ray photons into the model developed here. The present task was divided into two parts, namely, (1) the atomic scale model and (2) the macroscopic scale model, in which molecular dynamics (MD) and finite element method (FEM) simulation techniques were coupled. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers [20] (hereafter referred to as DA monomers), which were modelled using the MD technique. The degrees of SSP in the atomic scale upon different UV exposures were obtained in order to determine the absorption coefficients of the active layer, which were then input into the FEM model. The classical steady-state Helmholtz equation was engaged to model the re-emittance (i.e., reflection) from and transmittance through the active layer using the FEM technique. The multifrontal massively parallel sparse direct solver (MUMPS) [21] was employed to solve the present numerical problem.

Atomic Scale Model
To understand the degree of polymerization in an EBT3 film, we employed MD computations to model the behavior of DA monomers using the AIREBO (Adaptive Intermolecular Reactive Empirical Bond Order) potential [22]. The DA monomers were modeled in a rectangular simulation box with lengths of 800 Å in the x direction, 230 Å in the y direction, and 150 Å in the z direction, and the carboxylic chain (-COOH) was neglected. The number of oxygen atoms in the DA monomer was insignificant compared to those in the long hydrocarbon chain. The dominant bonding sites would be located along the simulated chain and the contribution of the -COOH group would be negligible. The attraction forces among the DA monomers were governed by attraction forces among atoms in the chains, which were inversely proportional to separations among the atoms. The simulation box consisted of ~10 6 (carbon + hydrogen) atoms in total, which accounted for ~15,000 DA monomer strands. A section of the simulation box is shown schematically in Figure 1.

Macroscopic Scale Model
The numerical model bundled with the MD code was based on the FEM technique, which simulated light propagation through an EBT3 film with different light-absorption coefficients. The FEM model simulated the realistic experimental set-up in which the EBT3 film was scanned. The classical steady-state Helmholtz equation was used to model the re-emittance (i.e., reflection) from and transmittance through the EBT3 film [23]: where φ was the intensity; S represented the source term; and D = 1/3(µ α + µ s ), µ α and µ s were absorption and reduced scattering coefficients, respectively. Equation (5) could be expanded as The Helmholtz's equation might be solved in 2-or 3-dimensional forms. Here it was solved in a 3-dimensional form that could directly emulate the experimental set-up. Equation (6) could be restated as the determination of function φ that minimized the functional χ defined as where Ω was the domain enclosing the volume. The Neumann condition was included in the functional, with the general form of the Neumann boundary condition as where n was the normal vector; and g and q represented the boundary flux and absorption term, respectively. The functional χ could be expressed as χ = χ f + χ s for convenience. Here, the FEM method was employed to solve the Helmholtz equation, which was based on dividing the region under analysis into a number of sub-regions (i.e., elements). The modelled domain (Ω) was a rectangular block with dimensions of 500 × 500 × 28 µm 3 . The subdivision of Ω was accomplished by introducing tetrahedron elements with 4 nodes (i.e., 4-noded elements). This technique is schematically summarized in Figure 2.
Appl. Sci. 2017, 7, 1031 4 of 13 (5) where ϕ was the intensity; S represented the source term; and D = 1/3(μα+μs), μα and μs were absorption and reduced scattering coefficients, respectively. Equation (5) could be expanded as (6) The Helmholtz's equation might be solved in 2-or 3-dimensional forms. Here it was solved in a 3-dimensional form that could directly emulate the experimental set-up. Equation (6) could be restated as the determination of function ϕ that minimized the functional χ defined as where Ω was the domain enclosing the volume. The Neumann condition was included in the functional, with the general form of the Neumann boundary condition as where n was the normal vector; and g and q represented the boundary flux and absorption term, respectively. The functional χ could be expressed as χ = χf + χs for convenience. Here, the FEM method was employed to solve the Helmholtz equation, which was based on dividing the region under analysis into a number of sub-regions (i.e., elements). The modelled domain (Ω) was a rectangular block with dimensions of 500 × 500 × 28 μm 3 . The subdivision of Ω was accomplished by introducing tetrahedron elements with 4 nodes (i.e., 4-noded elements). This technique is schematically summarized in Figure 2. Upon introduction of 4-noded elements, the dependent variable ϕ needed to be defined as a vector at each node as : The definition of ϕ would now be possible as (10) where i, j, l, m were nodes that emerged from the 4-nodes elements (see Figure 2). The variable ϕ could be expressed as , , , , , Upon introduction of 4-noded elements, the dependent variable φ needed to be defined as a vector at each node as The definition of φ would now be possible as where i, j, l, m were nodes that emerged from the 4-nodes elements (see Figure 2). The variable φ could be expressed as where N was the matrix containing the shape functions. For any set of values of {φ}, χ could be determined and minimization could be accomplished by writing the set of equations as Similarly, the elemental decomposition of χ would be considering χ e as the contribution of element e to the total integration. Interestingly, the relation between the quadratic function of the nodal values for each element would be linear: where the matrix S depended on the shape function N, which in fact was associated with the nodal elements as N = [N i , N j , N l , N m ]. The matrix S e was symmetric and could be represented as The term χ s was defined by Equation (16) could be redefined for other nodes of any element as Using the matrix form, it would be possible to state the total contribution for a 4-noded element as Equation (18) could be written as The two parts of the functional χ could now be combined as From Equations (12), (14), (19), The solution of Equation (21) gave the eigenvalues of µ α . This could be further extended for any number of associated nodes. Here, 4-noded tetrahedron elements were used. The mathematical description for a 4-noded tetrahedron element is given in Appendix A. The seamless mesh containing the tetrahedron elements is shown in Figure 3.
Appl. Sci. 2017, 7, 1031 6 of 13 (20) From Equations (12), (14), (19), The solution of Equation (21) gave the eigenvalues of μα. This could be further extended for any number of associated nodes. Here, 4-noded tetrahedron elements were used. The mathematical description for a 4-noded tetrahedron element is given in Appendix A. The seamless mesh containing the tetrahedron elements is shown in Figure 3. The employed mesh consisted of 468,345 domain elements including 34,056 boundary elements and 820 edge elements. The term D could be expressed in terms of absorption (μa) and reduced scattering (μs) coefficient as The present model was implemented and solved using the MUMPS solver, which allowed better exploitation of computer resources and reduced numerical fluctuations associated with the FEM technique. The solver configurations included: (1) a memory allocation factor of 1.2; (b) pivoting with a pivot threshold of 0.1; and (c) hyper-threading enabled for the multiprocessor cluster. Due to the employment of a fine mesh, the numerical simulations were executed on supercomputer with Intel ® Xeon E5-2670 2.60 GHz using 32 physical cores and hyper-threaded to 64.

Computation Scheme
The present approach to determine the net reflective optical density (Net ROD) (see below for definitions) is summarized in Figure 4. Variations in the absorption coefficient were obtained from the MD model, which would then be supplied to the FEM model. The intensity of reflection from the film surface would be scored. Ratios between different intensities were obtained to compare with experimentally determined Net ROD values. The employed mesh consisted of 468,345 domain elements including 34,056 boundary elements and 820 edge elements. The term D could be expressed in terms of absorption (µ a ) and reduced scattering (µ s ) coefficient as The present model was implemented and solved using the MUMPS solver, which allowed better exploitation of computer resources and reduced numerical fluctuations associated with the FEM technique. The solver configurations included: (1) a memory allocation factor of 1.2; (b) pivoting with a pivot threshold of 0.1; and (c) hyper-threading enabled for the multiprocessor cluster. Due to the employment of a fine mesh, the numerical simulations were executed on supercomputer with Intel ® Xeon E5-2670 2.60 GHz using 32 physical cores and hyper-threaded to 64.

Computation Scheme
The present approach to determine the net reflective optical density (Net ROD) (see below for definitions) is summarized in Figure 4. Variations in the absorption coefficient were obtained from the MD model, which would then be supplied to the FEM model. The intensity of reflection from the film surface would be scored. Ratios between different intensities were obtained to compare with experimentally determined Net ROD values.
Since Net ROD was a ratio of the reflected light intensity, the theoretical ratio Net ROD theoretical would not be affected by properties of the light source, including its strength, shape, size and distance from the film surface. For simplicity, a monoenergetic point source was chosen. The wavelength of light emission was assumed to be 650 nm. The irradiation of the film and the reflection from the film are schematically shown in Figure 5.

Computation Scheme
The present approach to determine the net reflective optical density (Net ROD) (see below for definitions) is summarized in Figure 4. Variations in the absorption coefficient were obtained from the MD model, which would then be supplied to the FEM model. The intensity of reflection from the film surface would be scored. Ratios between different intensities were obtained to compare with experimentally determined Net ROD values.  Since Net ROD was a ratio of the reflected light intensity, the theoretical ratio Net RODtheoretical would not be affected by properties of the light source, including its strength, shape, size and distance from the film surface. For simplicity, a monoenergetic point source was chosen. The wavelength of light emission was assumed to be 650 nm. The irradiation of the film and the reflection from the film are schematically shown in Figure 5.
where ϕunexposed and ϕexposed were reflected intensities of light from the unexposed and exposed film to UV radiation, respectively. The experimental Net ROD was defined as Net ROD = log(Pu/Pt), where Pu and Pt were the pixel values of the reflected intensities through an unexposed film and an exposed film, respectively. Here, the UV exposure values were taken from previous work by Chun et al. [8].

Results and Discussions
The absorption coefficients of EBT3 films versus UV exposure obtained from our model are shown in Figure 6. The increase in the absorption coefficient with increasing UV exposure could be explained in the atomic scale by the corresponding increase in the density of the system. The extra energy from larger UV exposures increased the kinetic energy of atoms in the monomer strands and facilitated their relocation to form a cross-linked structure. A denser polymerized configuration through solid-state polymerization would result from more bond formation. The lowest and highest absorption coefficients were found to be ~27.0 and ~29.0 cm −1 , respectively. The absorption coefficients were obtained from the product of mass energy-absorption coefficient with density. Here, the mass energy-absorption coefficient was assumed to be 100 cm 2 /g, but the chosen value was not important since Net ROD was the ratio of reflected intensities from the unexposed and exposed film to UV radiation. The variations in reflected light intensities from the EBT3 film were governed by the variation in density ρ versus UV exposure in other words, by the variation in absorption coefficient μα versus UV exposure, since μα = ρ × (μα)m. A snapshot from the atomic-scale system of monomer strands during relaxation after UV exposure and subsequent formation of bonds is shown in Figure 7. The Net ROD from the present model was calculated as where φ unexposed and φ exposed were reflected intensities of light from the unexposed and exposed film to UV radiation, respectively. The experimental Net ROD was defined as Net ROD = log(P u /P t ), where P u and P t were the pixel values of the reflected intensities through an unexposed film and an exposed film, respectively. Here, the UV exposure values were taken from previous work by Chun et al. [8].

Results and Discussions
The absorption coefficients of EBT3 films versus UV exposure obtained from our model are shown in Figure 6. The increase in the absorption coefficient with increasing UV exposure could be explained in the atomic scale by the corresponding increase in the density of the system. The extra energy from larger UV exposures increased the kinetic energy of atoms in the monomer strands and facilitated their relocation to form a cross-linked structure. A denser polymerized configuration through solid-state polymerization would result from more bond formation. The lowest and highest absorption coefficients were found to be~27.0 and~29.0 cm −1 , respectively. The absorption coefficients were obtained from the product of mass energy-absorption coefficient with density. Here, the mass energy-absorption coefficient was assumed to be 100 cm 2 /g, but the chosen value was not important since Net ROD was the ratio of reflected intensities from the unexposed and exposed film to UV radiation. The variations in reflected light intensities from the EBT3 film were governed by the variation in density ρ versus UV exposure in other words, by the variation in absorption coefficient µ α versus UV exposure, since µ α = ρ × (µ α ) m . A snapshot from the atomic-scale system of monomer strands during relaxation after UV exposure and subsequent formation of bonds is shown in Figure 7.  The reflected intensities (i.e., log(ϕ)) from the EBT3 films versus UV exposures are shown in Figure 8. A larger absorption coefficient led to increased absorption of light by the polymer material. For UV exposures larger than ~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was due to the insignificant variations in the absorption coefficient for UV exposures larger than ~40 J/cm 2 (see Figure 6). For relatively large UV exposures, the monomers became almost fully polymerized, and the densely cross-linked structure significantly lowered the efficiency for further bond formation between DA monomer strands. As a result, the absorption coefficient, which depended on the density of polymerization, also did not significantly increase at large UV exposures.  The reflected intensities (i.e., log(ϕ)) from the EBT3 films versus UV exposures are shown in Figure 8. A larger absorption coefficient led to increased absorption of light by the polymer material. For UV exposures larger than ~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was due to the insignificant variations in the absorption coefficient for UV exposures larger than ~40 J/cm 2 (see Figure 6). For relatively large UV exposures, the monomers became almost fully polymerized, and the densely cross-linked structure significantly lowered the efficiency for further bond formation between DA monomer strands. As a result, the absorption coefficient, which depended on the density of polymerization, also did not significantly increase at large UV exposures. The reflected intensities (i.e., log(φ)) from the EBT3 films versus UV exposures are shown in Figure 8. A larger absorption coefficient led to increased absorption of light by the polymer material. For UV exposures larger than~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was due to the insignificant variations in the absorption coefficient for UV exposures larger than~40 J/cm 2 (see Figure 6). For relatively large UV exposures, the monomers became almost fully polymerized, and the densely cross-linked structure significantly lowered the efficiency for further bond formation between DA monomer strands. As a result, the absorption coefficient, which depended on the density of polymerization, also did not significantly increase at large UV exposures. In the model here, Net ROD was defined as Net ROD = log(ϕunexposed/ϕexposed), where ϕunexposed and ϕexposed were reflected intensities of light from the unexposed and exposed film to UV radiation, respectively, in line with the definition of experimental Net ROD as Net ROD = log(Pu/Pt), where Pu and Pt were the pixel values of reflected intensities through an unexposed film and an exposed film, respectively [8]. It is remarked here that, during the conversion of experimentally obtained images of the EBT3 film into grayscale, the pixel values ranged from 255 (completely white) to 0 (completely black), so the ratio Pu/Pt was always larger than unity. Similarly, the ratio ϕunexposed/ϕexposed ratio was always larger than unity. Both pixel values and reflected light intensities could be traced back to properties of the EBT3 film under analysis.
The comparisons among Net ROD values (obtained experimentally by Chun et al. [8]) and from theoretical predictions with two different mass energy-absorption coefficients, (μα)m versus UV-light exposures are shown in Figure 9. The Net ROD values from theoretical predictions were scored over the entire film surface that had been exposed to UV. In the FEM technique, an element consisted of a number of nodes (4 nodes for the tetrahedron used in the present work) each being associated with a shape function. The intensity of the reflected light in an element was determined through summing the intensities over all nodes with consideration of their corresponding shape function. The average reflected intensity was determined as the arithmetic mean of the intensities of elements in the modeled domain [24]. The theoretically obtained Net ROD values were generally in good agreement with the experimental data and, as described above, did not vary noticeably with the chosen mass energy-absorption coefficients (viz., 100 vs 10 cm 2 /g). The two curves were close to each other, but the numerical values were not exactly the same. When different mass energy-absorption coefficients were employed, the reflected light intensities from the unexposed and exposed films changed by almost the same factor, so the Net ROD values were close. Figure 9 shows that Net ROD increases with the UV exposure. Upon UV exposures, the DA monomer strands gained enough energy to form a dense cross-linked polymerized structure, so the EBT3 film became darker. However, upon excessive UV exposures, polymerization would become saturated and there would be no further film darkening.
The deviations between the model and the experimental data in Figure 9 could be due to assumptions and simplifications made in the present model; such as, no other elements than carbon and hydrogen atoms were present in the EBT3 film, and no light leaked from the film. The spatial distribution of Net ROD over the upper and lower surface of an EBT3 film after being exposed to 164.70 J/cm 2 UV radiation is shown in Figure 10. The Net ROD values were larger at the upper surface, which was explained by the loss of intensities of light rays during their traversals of the active layer of the EBT3 film (thickness = 28 μm). The largest theoretical Net ROD was recorded at the center of the film because a point source was used in the simulations. In the model here, Net ROD was defined as Net ROD = log(φ unexposed /φ exposed ), where φ unexposed and φ exposed were reflected intensities of light from the unexposed and exposed film to UV radiation, respectively, in line with the definition of experimental Net ROD as Net ROD = log(P u /P t ), where P u and P t were the pixel values of reflected intensities through an unexposed film and an exposed film, respectively [8]. It is remarked here that, during the conversion of experimentally obtained images of the EBT3 film into grayscale, the pixel values ranged from 255 (completely white) to 0 (completely black), so the ratio P u /P t was always larger than unity. Similarly, the ratio φ unexposed /φ exposed ratio was always larger than unity. Both pixel values and reflected light intensities could be traced back to properties of the EBT3 film under analysis.
The comparisons among Net ROD values (obtained experimentally by Chun et al. [8]) and from theoretical predictions with two different mass energy-absorption coefficients, (µ α ) m versus UV-light exposures are shown in Figure 9. The Net ROD values from theoretical predictions were scored over the entire film surface that had been exposed to UV. In the FEM technique, an element consisted of a number of nodes (4 nodes for the tetrahedron used in the present work) each being associated with a shape function. The intensity of the reflected light in an element was determined through summing the intensities over all nodes with consideration of their corresponding shape function. The average reflected intensity was determined as the arithmetic mean of the intensities of elements in the modeled domain [24]. The theoretically obtained Net ROD values were generally in good agreement with the experimental data and, as described above, did not vary noticeably with the chosen mass energy-absorption coefficients (viz., 100 vs. 10 cm 2 /g). The two curves were close to each other, but the numerical values were not exactly the same. When different mass energy-absorption coefficients were employed, the reflected light intensities from the unexposed and exposed films changed by almost the same factor, so the Net ROD values were close. Figure 9 shows that Net ROD increases with the UV exposure. Upon UV exposures, the DA monomer strands gained enough energy to form a dense cross-linked polymerized structure, so the EBT3 film became darker. However, upon excessive UV exposures, polymerization would become saturated and there would be no further film darkening.
The deviations between the model and the experimental data in Figure 9 could be due to assumptions and simplifications made in the present model; such as, no other elements than carbon and hydrogen atoms were present in the EBT3 film, and no light leaked from the film. The spatial distribution of Net ROD over the upper and lower surface of an EBT3 film after being exposed to 164.70 J/cm 2 UV radiation is shown in Figure 10. The Net ROD values were larger at the upper surface, which was explained by the loss of intensities of light rays during their traversals of the active layer of the EBT3 film (thickness = 28 µm). The largest theoretical Net ROD was recorded at the center of the film because a point source was used in the simulations. (a) (b) Figure 10. Spatial distribution of Net ROD over (a) upper and (b) lower surface of EBT3 film after being exposed to 164.70 J/cm 2 UV radiation.

Conclusions
In the present work, we presented the first numerical model for the simulation of coloration changes in radiochromic EBT3 film. A MD-coupled FEM model was developed to study the coloration change of the EBT3 film when exposed to UV radiation. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using MD. The densities of the active layer of the EBT3 film were determined upon different UV exposures. The results obtained from MD computations were used as inputs to the FEM model to simulate light propagation through the EBT3 film. The classical steady-state Helmholtz equation was implemented in the FEM model to determine the reflected light intensities. Very good agreement between experimentally and theoretically obtained coloration in terms of net reflective optical density was achieved for different UV exposures. In particular, for UV exposures larger than ~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization and, thus, lower efficiency of further bond formation between DA monomer strands. With the success of our present model for explaining solid-state polymerization (SSP) in the EBT3 film, we are in a good position to perform related studies on dose responses for X-ray photons by adding processes involving secondary electrons created by X-ray photons, which will be covered in a future paper. (a) (b) Figure 10. Spatial distribution of Net ROD over (a) upper and (b) lower surface of EBT3 film after being exposed to 164.70 J/cm 2 UV radiation.

Conclusions
In the present work, we presented the first numerical model for the simulation of coloration changes in radiochromic EBT3 film. A MD-coupled FEM model was developed to study the coloration change of the EBT3 film when exposed to UV radiation. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using MD. The densities of the active layer of the EBT3 film were determined upon different UV exposures. The results obtained from MD computations were used as inputs to the FEM model to simulate light propagation through the EBT3 film. The classical steady-state Helmholtz equation was implemented in the FEM model to determine the reflected light intensities. Very good agreement between experimentally and theoretically obtained coloration in terms of net reflective optical density was achieved for different UV exposures. In particular, for UV exposures larger than ~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization and, thus, lower efficiency of further bond formation between DA monomer strands. With the success of our present model for explaining solid-state polymerization (SSP) in the EBT3 film, we are in a good position to perform related studies on dose responses for X-ray photons by adding processes involving secondary electrons created by X-ray photons, which will be covered in a future paper. Figure 10. Spatial distribution of Net ROD over (a) upper and (b) lower surface of EBT3 film after being exposed to 164.70 J/cm 2 UV radiation.

Conclusions
In the present work, we presented the first numerical model for the simulation of coloration changes in radiochromic EBT3 film. A MD-coupled FEM model was developed to study the coloration change of the EBT3 film when exposed to UV radiation. The active layer in the EBT3 film consisted of diacetylene (DA) pentacosa-10,12-diynoate monomers, which were modelled using MD. The densities of the active layer of the EBT3 film were determined upon different UV exposures. The results obtained from MD computations were used as inputs to the FEM model to simulate light propagation through the EBT3 film. The classical steady-state Helmholtz equation was implemented in the FEM model to determine the reflected light intensities. Very good agreement between experimentally and theoretically obtained coloration in terms of net reflective optical density was achieved for different UV exposures. In particular, for UV exposures larger than~40 J/cm 2 , the reflected light intensity decreased at a lower rate when compared to other UV exposure values, which was explained by the densely cross-linked structure under near-complete polymerization and, thus, lower efficiency of further bond formation between DA monomer strands. With the success of our present model for explaining solid-state polymerization (SSP) in the EBT3 film, we are in a good position to perform related studies on dose responses for X-ray photons by adding processes involving secondary electrons created by X-ray photons, which will be covered in a future paper.
Tetrahedron element: The shape function of the tetrahedron element can be expressed in the form of a linear polynomial as For the 4 different nodes, the shape function can be decomposed into nodal components as where i, j, l and m are the nodes from the 4-noded tetrahedron element. Therefore, φ can be expressed as a m , however further derivation of these terms is not shown in the present paper. A distinctive term in the matrix S can now be represented as where V is the volume of the 4-noded tetrahedron element.