Mathematical optimization approach for estimating the quantum yield distribution of a photochromic reaction in a polymer

The convolution of a series of events is often observed for a variety of phenomena such as the oscillation of a string. A photochemical reaction of a molecule is characterized by a time constant, but materials in the real world contain several molecules with different time constants. Therefore, the kinetics of photochemical reactions of the materials are usually observed with a complexity comparable with those of theoretical kinetic equations. Analysis of the components of the kinetics is quite important for the development of advanced materials. However, with a limited number of exceptions, deconvolution of the observed kinetics has not yet been mathematically solved. In this study, we propose a mathematical optimization approach for estimating the quantum yield distribution of a photochromic reaction in a polymer. In the proposed approach, time-series data of absorbances are acquired and an estimate of the quantum yield distribution is obtained. To estimate the distribution, we solve a mathematical optimization problem to minimize the difference between the input data and a model. This optimization problem involves a differential equation constrained on a functional space as the variable lies in the space of probability distribution functions and the constraints arise from reaction rate equations. This problem can be reformulated as a convex quadratic optimization problem and can be efﬁciently solved by discretization. Numerical results are also reported here,

The convolution of a series of events is often observed for a variety of phenomena such as the oscillation of a string. A photochemical reaction of a molecule is characterized by a time constant, but materials in the real world contain several molecules with different time constants. Therefore, the kinetics of photochemical reactions of the materials are usually observed with a complexity comparable with those of theoretical kinetic equations. Analysis of the components of the kinetics is quite important for the development of advanced materials. However, with a limited number of exceptions, deconvolution of the observed kinetics has not yet been mathematically solved. In this study, we propose a mathematical optimization approach for estimating the quantum yield distribution of a photochromic reaction in a polymer. In the proposed approach, time-series data of absorbances are acquired and an estimate of the quantum yield distribution is obtained. To estimate the distribution, we solve a mathematical optimization problem to minimize the difference between the input data and a model. This optimization problem involves a differential equation constrained on a functional space as the variable lies in the space of probability distribution functions and the constraints arise from reaction rate equations. This problem can be reformulated as a convex quadratic optimization problem and can be efficiently solved by discretization. Numerical results are also reported here, and they verify the effectiveness of our approach. © 2017 Author(s). All

I. INTRODUCTION
Photochromism is a chemical reaction between two isomers of a molecule. It is induced by light absorption, which reversibly converts one isomer to another. 1 One example of a photochromic reaction is the trans to cis isomerization of azobenzene, as shown in Figure 1.
A structural change in the molecule is accompanied by changes in various physical properties such as color, refractive index, polarity, wettability, and electric conductivity. These changes in properties are exploited as functional attributes of the materials and are induced by photo-irradiation as an external switch. Examples of the applications of photochromic molecules are photon-mode optical data storage, color-modulatable glasses, photo-optical materials, photo-mechanical materials, and self-organized structured solar-cell materials. 2 One of the physical constants of photochromic reactions is quantum yield, which is defined as the ratio of the number of reacted molecules to the number of absorbed photons, and which is usually unique to a molecule. Most organic photo-functional materials comprise photochromic molecules dispersed in a solid-state polymer matrix, where there is a distribution of free volume or the size of the space between polymer chains. Hence, the environment of individual photochromic molecules is different from each other. Therefore, the quantum yield of a photochromic reaction in a solid-state polymer has a large distribution that reflects the inhomogeneity of the polymer matrix, even though the quantum yield of the reaction in a solution is a unique value without a distribution.
There have been several attempts to analyze the multicomponent kinetics. A trivial example is Fourier transformation analysis. One can precisely determine the distribution of the components for a given wave if the observed event is known as the convolution of sine functions. However, most natural events cannot be mathematically solved because their kinetics are mutually coupled with no mutual orthogonal basis set for each.
Representative kinetics for a single molecular reaction is expressed as an exponentially decaying function characterized by a time constant, following the separation of variables that describe the differential equation. To express the inhomogeneity of the time constant, a reaction in viscous media may be analyzed with an extended exponential equation introducing another parameter in addition to the original time constant. 8 Another approach to the analysis of multicomponent kinetics is the deconvolution method, which is usually performed for fluorescence lifetime measurements. 9 One can determine the best-fit solution for a set of parameters by assuming that the fluorescence intensity I can be expressed as where τ 1 and τ 2 are the lifetimes of governing the decaying components of fluorescence; A 1 and A 2 are the pre-exponential factors, respectively. The problem is that only a limited number of components of the reaction can be expressed on the bases of past methodologies even though there should be a wide distribution of components with an unknown distribution function in practical materials. We analyzed for the first time the quantum yield distribution of the photochromic reaction of azobenzene in a solid-state polymer to determine several valuable principles of molecular design for high performance materials. However, to simplify the calculations, the analysis was performed assuming that the quantum yield follows a Gaussian distribution. Further progress in the analysis is still limited. 10,11 Therefore, the development of a novel methodology for solving the distribution function of the components for given data of the convolution of complex phenomena without any assumptions is very important in material science.
In this study, we propose a technique for obtaining the quantum yield distribution of a photochromic reaction based on the time-dependent change in absorbance of a sample. This change in absorbance is described in theory as a convex combination of decaying curves with their respective quantum yields. However, the experimental data also contains noise arising for example from the fluctuation in light intensity, electric noise from detectors, and light scattering at the surface of sample films. Therefore, the experimental data cannot be fitted solely by the convolution of theoretical equations because of noise, which makes it difficult to estimate the true distribution of the theoryderived equation in the experimental data. Because of this, we propose to solve an optimization problem to obtain a density function that minimizes the difference between the measured and theoretical absorbances corresponding to the density function and to regard it as an estimate. However, it is difficult to solve the optimization problem directly because it is an optimization problem over a functional space and has a differential equation as constraint. In contrast, it is sufficient to obtain a discretized distribution instead of a continuous distribution. Therefore, we discretized the quantum yield distribution to overcome these difficulties. From the discretization, the estimation of the quantum yield density function is reduced to the computation of the quantum yield of each sample point. In addition, we can remove the differential equation constraints by solving the differential equation for each sample point, without loss of equivalence. The resulting optimization model is a convex quadratic optimization problem and is thus easily solved with an existing solver. An optimal solution to the discretized problem is close to that of the original problem if we discretize it sufficiently finely. In general, the solution of the discretized problem converges to the solution of the original problem as mentioned in Still. 12 The remainder of this paper is organized as follows: In Section II, we derive a reaction rate equation for photochromism. In Section III, we cover the proposed optimization model. In Section IV, we report on the results of numerical results. Section V is devoted to concluding remarks.

II. REACTION RATE FOR A UNIFORM QUANTUM YIELD
We briefly mention the reaction kinetics of photochromism. First, we introduce the Beer-Lambert law. Then, we derive a reaction rate equation for photochromic molecules. More details about the reaction kinetics of photochromism can be found in the physical chemistry textbook. 13 The Beer-Lambert law is a formula of the absorption of light by a material through which light is traveling. Consider light with intensity I 0 [µmol · m −2 · s −1 ] irradiating a polymer sample of optical path length l [cm]. It is known that the intensity I of the transmitted light decreases exponentially with increasing l. That is, I = I 0 e al or equivalently where the constant b = a log 10 e is called the coefficient of absorption. The left-hand side log 10 (I 0 /I) is called the absorbance. In addition, we assume there is no interaction between photochromic molecules. Then, the coefficient absorption b is proportional to the molar concentration c [mol · L −1 ] of photochromic molecules. Let the constant of proportionality be [mol −1 · L · cm −1 ]. That is called the molar coefficient of absorption. Then, Equation (1) is reformulated as which is referred to as the Beer-Lambert law. Let us consider a situation where the the molar concentration changes sdepending with time. Let c(t) be the molar concentration at time t. Then, the absorbance at time t is denoted by X(t) = c(t)l. From the Beer-Lambert law, we see that the number of photons absorbed by the sample at time t is (I 0 I)S = I 0 S(1 10 X (t ) ), where S is the cross-sectional area of the sample. Next, we derive a reaction rate equation for photochromic molecules. Here we consider transisomers and cis-isomers photochemically changing from one to the other, such as with the photochromism of azobenzene. In this situation, the absorbance X(t) becomes the sum of that of the trans-isomers and the cis-isomers. Hence, we obtain where trans and cis denote the molar coefficients of absorption for the trans-isomer and the cis-isomer, respectively; similarly, c trans (t) and c cis (t) denote the molar concentrations of the trans-isomer and the cis-isomer at time t, respectively. From the results of the Beer-Lambert law, the number of absorbed photons for the trans-isomers is written as trans c trans (t)l cis c cis (t)l + trans c trans (t)l I 0 S(1 − 10 −X(t) ).
In addition, let φ trans→cis be the quantum yield of the reaction trans-isomers to cis-isomers and φ cis→trans be that of the reverse reaction. Then, we can obtain the transformation rate from trans-isomer to cis in [mol · s −1 ] as follows d(10 −3 c trans (t)Sl) dt = − trans c trans (t)l cis c cis (t)l + trans c trans (t)l cis c cis (t)l cis c cis (t)l + trans c trans (t)l The left-hand side is the definition of the transformation rate. The first term on the right-hand side is the decreasing rate of the trans-isomers for the forward reaction and the second is the increasing rate of the trans-isomers for the backward reaction. In our setting, c cis (0) = 0 and hence for any t ≥ 0 hold. Using basic calculus (see Appendix), we obtain the following ordinary differential equation (ODE) of the absorbance X(t) where r φ = φ cis→trans /φ trans→cis and r = cis / trans are substance-specific constants. A solution to Equation (5) depends on the value of the quantum yield φ trans→cis of the forward reaction. In what follows, let X(t; φ) be a solution to Equation (5) with φ trans→cis = φ. Hence, if function X(t; φ) satisfies Equation (5) with φ trans→cis = φ, we denote the solution as F(t, X(t; φ),Ẋ(t; φ); φ) = 0. We plot the numerical solution to Equation (5) using various values of the quantum yield φ and the measured absorbance (y(t) : t ∈ T ) in Figure 2, where T is the set of measured times. Details regarding the sample are covered in Section IV. The measured data are plotted as a spiky line. The spikes in the measured data can be regarded as deriving from noise. In contrast, the numerical solutions are drawn as smooth curves. They monotonically decrease and converge to a constant, except at φ = 0, the horizontal dotted line. The convergence corresponds to the equilibrium of the forward and backward reactions. Note that the curve of the measured data (y(t) : t ∈ T ) does not match any curve from the numerical solution X(t; φ). This implies that the quantum yield of the sample is not constant but has a distribution instead.

III. OPTIMIZATION MODEL FOR ESTIMATING QUANTUM YIELD DISTRIBUTION
We now consider an optimization model to estimate the density function of the quantum yield distribution from the data of the measured absorbance and the reaction rate in Equation (5). Specifically, we propose a method for finding the most suitable density function β for the given measured data ( y(t) : t ∈ T ) and the parameters in Equation (5). Here, β(φ) denotes the probability density of the existence of molecules whose quantum yield φ trans→cis of the forward reaction is φ. In what follows, we assume that 0 ≤ φ trans→cis ≤ 1. It is equivalent to assuming that no chain reaction occurs during photochromism.
If the quantum yield of the forward reaction takes a significant value φ, we can assume that the measured absorbance y(t) at each time t ∈ T is a solution X(t; φ) of Equation (5) . Similarly, if the quantum yield has a distribution, we can naturally assume that the measured absorbance y(t) at time t ∈ T is an expectation of the solution X(t; φ) to Equation (5) with the density function β and noise ε t . That is, This model assumes that the measured data (y(t) : t ∈ T ) is a noisy convex combination of X(t; φ)'s with weight β(φ).
Under this assumption, let us consider a way to obtain the most suitable probability distribution function β for the measured data. We can model the problem to minimize the sum of squares of the noise intensities (ε t : t ∈ T ) under X(t; φ) to satisfy an initial state condition and Equation (5) as follows: This optimization model can be regarded as a constrained least-squares method or a constrained linear regression model by regarding (y(t) : t ∈ T ) as explained variables, (X(t; φ) : t ∈ T ) as explanatory variables, and β(φ) (0 ≤ φ ≤ 1) as coefficients. The optimization problem expressed by Equation (6) is ill-conditioned. Indeed, for two efficiently close quantum yields, say φ and φ , X(t; φ) is close to X(t; φ ). Hence, there are multiple feasible solutions that attain the same objective value. In the discretized model, to be described below, we see that the matrix of coefficients for the quadratic term in the objective function is ill-conditioned. In terms of linear regression, this is referred to as multicollinearity.
We employ the so-called regularization technique to overcome this problem. In our situation, the L 1 -regularization is not suitable because we cannot invoke the sparsity in elaborating a solution. Here, we employ the standard L 2 -regularization. Let λ (≥ 0) be a parameter to define the strength of the regularization. By introducing the L 2 -regularization term to the objective function of the optimization problem expressed by Equation (6), we obtain the following optimization problem: This optimization model may be regarded as a constrained ridge regression. The optimization problem expressed by Equation (7) is an ODE-constrained optimization problem on a functional space. Hence, this problem is difficult to solve. Hence, in what follows, by discretizing the problem, we reformulate the model into a convex quadratic programming (QP) model. In more concrete terms, using sampling points φ 0 < · · · < φ N , we approximate the integral There are multiple means to discrete an integral. Using either Gauss quadrature or Clenshaw-Curtis quadrature, 14,15 etc. we can achieve an approximation with high-accuracy. In addition, by numerically solving F(t, X(t; φ i ),Ẋ(t; φ i )) = 0 with the initial condition X(0; φ i ) = X 0 for i = 0, . . ., N, we can compute a numerical solution (X(t; φ i ) : t ∈ T ). By substituting the numerical solution into the optimization problem expressed by Equation (7), we obtain the following optimization model: As this model is only a convex QP problem, we can easily solve a moderate-sized instance by employing off-the-shelf solvers. [16][17][18]

IV. NUMERICAL RESULTS
We now report our numerical results that verify the effectiveness of our approach. Our numerical experiments comprise two parts. One was an experiment to estimate the quantum yield distribution using artificial data generated from a given true quantum yield distribution. We verified that the proposed method estimates with high accuracy the true distribution. The other used empirical data obtained from our measurements. We determined the unknown quantum yield distribution of the sample used in the experiments.
We executed all our experiments in the following environment:

A. Numerical results for artificial data
We performed the experiments using the following procedure: 1. Choose a true quantum yield distribution β(φ).
2. Choose sampling points φ i for i = 1, . . ., N and numerically solve the reaction rate equation F(t, X(t; φ),Ẋ(t; φ i ); φ i ) = 0 with an appropriate initial condition. 3. Define the set T of (virtual) measuring times, numerically integrate y(t) = ∫ 1 0 X(t; φ) β(φ)dφ for all t ∈ T , and treat them as artificial data of absorbance. 4. For multiple values of the regularization parameter λ, solve the L 2 -regularized model, i,e., Equation (8), and compute the difference between the optimal solutionβ and the true distribution β.
We used the following distributions as true quantum yield distributions: Gaussian distributions N(0.05, 0.01 2 ) and N(0.05, 0.005 2 ), beta distributions Be (1,19) and and Mix B = 0.6N(0.03, 0.005 2 ) + 0.4N(0.08, 0.005 2 ). We set the set of measured times to T = {0, 30, 60, . . ., 55320} the same as that for the empirical data mentioned in Section IV B. We used the trapezoidal rule in the numerical integration and discretization of integrals. We measured the difference between the optimal solutionβ and the true distribution β using numerical integration of the Kullback-Leibler divergence 19 of β forβ. That is, We solved instances of the optimization problem expressed by Equation (8) with multiple values of the regularization parameter λ. The change in the difference betweenβ and β is shown in Figure 3. We see the following trend for all true distributions: when the regularization parameter λ is increased, the difference D( β||β) between the optimal solutionβ and the true distribution β decreases until reaching a minimum near λ = 10 −9 ; thereafter, it increases.
Next, for each true distribution, we plotted the true distribution (true), an optimal value with λ = 0 (naive), and an optimal value with an appropriate λ that minimized D( β||β) (appropriate) in  from the true distribution. In contrast, the optimal solution with an appropriate λ is close to the true distribution. Indeed, the two lines overlap in Figure 6. These results imply the following: The regularization parameter λ has an appropriate value; By employing our proposed approach with an appropriate λ, we can estimate with high accuracy the quantum yield distribution.

B. Numerical results from empirical data
In the numerical experiments with empirical data, we used the following data measured with the following procedures. Sample films were prepared by bulk polymerization of methyl methacrylate containing azobenzene using AIBN as an initiator. The change in the absorbance of the sample film was recorded by measuring the transmitted light intensity with a photodiode under photoirradiation with 365 nm light from a Xe lamp using a monochromator every 30 s until 55320 s elapsed. The quantum yield of this reaction is known not to exceed unity. We show the optimal solutions in Figure 10 for when we used the regularization parameters λ = 10 −9 , 10 −7 , 10 −5 . These values were determined from the numerical results in Section IV A. For each case, we obtained a distribution that has two peaks near φ = 0 and φ = 0.055. The peak at φ = 0 implies the existence of molecules that cannot react because of lack of free volume. In contrast, if a molecule has sufficient free volume, photochromism occurs with quantum yield φ = 0.055.

V. CONCLUDING REMARKS
We proposed an optimization model to estimate the quantum yield distribution of photochromic molecules in a polymer. The problem of finding the most suitable probability distribution function for the quantum yield distribution was modeled as a differential equation constrained optimization problem on a functional space. This can be reduced to a convex quadratic optimization problem by discretizing the space of the quantum yield. Our numerical results imply that we can estimate with high accuracy the quantum yield distribution by solving the proposed optimization problem with an appropriate regularization parameter, and we can quantitatively analyze empirical data.
In our numerical experiments, we defined an appropriate regularization parameter by executing preparatory experiments. A more objective way is yet to be determined. In terms of material science, our future work includes the development of new devices. To achieve this objective, software based on our approach must be implemented. In addition, a way to solve the optimization problem expressed by Equation (7) without discretization and one to calculate interval estimate of the distribution is an open problem for investigation.