Role of heat generation and thermal diffusion during frontal photopolymerization

Frontal photopolymerization (FPP) is a rapid and versatile solidification process that can be used to fabricate complex three-dimensional structures by selectively exposing a photosensitive monomer-rich bath to light. A characteristic feature of FPP is the appearance of a sharp polymerization front that propagates into the bath as a planar traveling wave. In this paper, we introduce a theoretical model to determine how heat generation during photopolymerization influences the kinetics of wave propagation as well as the monomer-to-polymer conversion profile, both of which are relevant for FPP applications and experimentally measurable. When thermal diffusion is sufficiently fast relative to the rate of polymerization, the system evolves as if it were isothermal. However, when thermal diffusion is slow, a thermal wavefront develops and propagates at the same rate as the polymerization front. This leads to an accumulation of heat behind the polymerization front which can result in a significant sharpening of the conversion profile and acceleration of the growth of the solid. Our results also suggest that a novel way to tailor the dynamics of FPP is by imposing a temperature gradient along the growth direction.


I. INTRODUCTION
Photopolymerization is a robust solidification process that takes place when a photosensitive monomer-rich bath is exposed to light. The absorption of radiation by this bath initiates a series of chemical reactions that ultimately leads to polymerization, crosslinking of a polymer network, or both. The high rates of polymerization that are obtained by exposing the bath to high-intensity radiation facilitate the rapid fabrication of complex three-dimensional network solids when multifunctional monomers are employed. The pattern of the resulting network can be precisely controlled through selective illumination of the bath and the development (or selective removal) of unpolymerized material. The lightsensitive monomer thus behaves as a "negative photoresist" in the context of lithographic processes. Photopolymerization is particularly attractive from a manufacturing perspective because it can be performed at room temperature and in open or controlled atmospheres using readily available UV light sources and commercial monomer or prepolymer formulations [1][2][3]. In practice, photopolymerization has been used in a number of contexts [3], ranging from photolithography [4,5], three-dimensional stereolithography [6], rapid prototyping [7,8], coatings [9], adhesives [10], biomedicine [11], tissue engineering [12,13], and dentistry [14], all of which demonstrate the great versatility of this fabrication method as well as the exciting possibility of manufacturing with light.
The absorption of incoming radiation by the monomer-rich bath leads to the intensity decaying with distance from the illuminated surface, as expected from the Beer-Lambert law [ Fig. 1(a)]. The gradient in radiation intensity has important consequences for the solidification process, as it drives the directional growth of the network. Due to the overall rate of polymerization being proportional to the intensity, the network solid grows from the illuminated surface towards the bath in a direction that follows the intensity gradient [ Fig. 1(b)]. That is, the network grows from regions of high intensity to low intensity. Under specific experimental conditions that will be discussed below, the interfacial layer separating solid-rich and liquid-rich phases can be thin in comparison to the typical size of the solid. In this case, the rapid transition from solid to liquid is analogous to a sharp solidification front and characteristic of a frontal photopolymerization (FPP) process [4,15,16]. FPP is, thus, defined by the presence of a sharp polymerization front that propagates into the bath as a traveling wave, leading to a scenario that is qualitatively similar to traditional directional solidification [17].
The evolution of the front position can be modelled by first introducing an order parameter φ measuring the local fraction of monomer that has been converted into polymer [ Fig. 1(c)]. The polymerization front, corresponding to the leading edge of the solid, can be formally defined as the location where the order parameter reaches a critical value φ c . The position of the front is, therefore, implicitly defined through φ rather than determined from an explicit equation of motion.
Practical implementations of FPP greatly benefit from quantitative knowledge of how the shape of the monomer-topolymer conversion profile, as well as the propagation of the polymerization front, depend on the experimental parameters. By determining the relationship between exposure time t and the position of the polymerization front, i.e., the height of the solid [ Fig. 1(d)], Cabral et al. [4] demonstrated that complex, multilevel patterns (used as microfluidic devices) can be rapidly produced with FPP. The material properties of the cured solid, e.g., its elastic modulus [18] or refractive index, can be strongly dependent on the conversion fraction φ; therefore, the ability to finely tune the conversion profile is directly relevant to the fabrication of gradient polymer solids having material properties that vary in a systematic manner. Indeed, as demonstrated by Turturro et al. [12,13], FPP offers a flexible approach for the controlled production of hydrogel scaffolds with gradient material properties used in tissue engineering. The manufacturing potential of FPP has stimulated a number of works focusing on the development of theoretical models [16,[19][20][21][22][23] that can be used as predictive tools, as well as experimental research [15,18,24,25] that aims to identify the key mechanisms underlying photopolymerization. In previous studies [4,15,18], we have proposed that FPP requires (i) strong optical attenuation with limited (ii) mass and (iii) thermal transport. Recently [26], we have confirmed and established quantitatively that conditions (i) and (ii) are, indeed, necessary for FPP to take place. One of the purposes of this paper is to determine whether condition (iii) is also a firm requirement by systematically investigating, from a theoretical point of view, the relationship between heat generation during photopolymerization and the onset and dynamics of FPP. We will pay particular attention to understanding how the motion of the polymerization front and the shape of the conversion profile are influenced by thermal effects, as these quantities are both relevant for FPP applications and experimentally measurable [4,18].
In many photopolymerizing systems, exothermic reactions are a significant source of heat generation [27,28], with the continual absorption of radiation by the medium also being a contributing factor [29]. The Arrhenius-like dependence of various reaction rates introduces a coupling between heat generation and polymerization, enabling these two mechanisms to interact with each other in a nonlinear fashion. The main purpose of this paper is, therefore, to investigate how thermal effects impact front propagation and determine their consequences for practical manufacturing processes.
It is important to emphasize that despite the coupling between photopolymerization and temperature, FPP remains distinct from thermal frontal polymerization (TFP). Whereas TFP is driven solely by the positive feedback loop arising from the release of thermal energy, the onset of thermal polymerization, and Arrhenius kinetics, FPP is driven by the absorption of radiation and can be terminated by stopping the illumination. A related process is isothermal frontal polymerization (IFT). Like TFP, IFP is autocatalytic; however, IFP is a consequence of the Tromsdorff-Norish effect, whereby polymerization inhibits terminating reactions by increasing the local viscosity of the mixture. For a comprehensive review of these frontal processes, the reader is referred to Pojman [30].
Recent interest in photopolymerization has sparked the development of a range of models that can be roughly grouped into two categories depending on their complexity. In one group are physicochemical models that accurately account for each of the reaction steps, e.g., photolysis, photoinitiation, propagation, chain transfer, and termination [31][32][33][34]; oxygen inhibition [35,36]; nonuniform distributions in the length of polymer chains [37,38]; heat generation and transport [20]; mass transport due to convection and/or diffusion [25,39]; and optical effects such as scattering [40] and refractive index modulation [41]. Although such models can offer theoretical insights into photopolymerization processes, their practical use is limited by a lack of tractability and large number of parameters, some of which cannot be measured experimentally.
Overcoming these drawbacks has motivated the development of a second group of models, the "minimal" or coarse-grained models. These make use of phenomenological, rather than explicit, physicochemical descriptions of photopolymerization and are specifically designed so they only depend on, and capture the evolution of, experimental observables, e.g., the height of the growing solid and radiation intensity. In light of this, we opt to extend the minimal model introduced by Cabral et al. [4], which is able to accurately capture the experimental observations of a range of (near isothermal) thiol-ene systems with low values of φ c [4,15,18,26].
In Sec. II of this paper, we outline some key results from the minimal model of Cabral et al. [4] describing isothermal FPP. These will serve as baseline solutions that will be compared to those from an extended thermal model accounting for heat generation and transport. In Sec. III we present the thermal model and a detailed analysis. The results are discussed in a practical context and the paper concludes in Sec. IV.

II. PROBLEM DEFINITION AND OVERVIEW OF PREVIOUS RESULTS
We consider the physical situation shown in Fig. 1, whereby a photosensitive monomer-rich bath is placed on a substrate with a prescribed temperature and covered by a transparent surface, e.g., a glass slide. The bath is illuminated by a light source with a given intensity, thus driving photopolymerization and the growth of a network solid from the transparent surface into the bulk.
A simple model to describe this process has been proposed by Cabral et al. [4], which considers the isothermal evolution of an order parameter φ = φ(z,t) representing the local fraction of polymer at a point z in the mixture at time t, and the intensity of radiation throughout the mixture, I = I (z,t). The coordinate system is chosen such that z points downwards with z = 0 corresponding to the bottom of the illuminated surface; see Fig. 1. We will often refer to φ as the conversion fraction and the curve φ = φ(z,t * ) for a fixed t * as the conversion profile. The rate of conversion is assumed to be proportional to an effective, overall, rate constant K with units of m 2 /J, the fraction of available monomer (1 − φ), and the intensity of radiation, I . Mass transport mechanisms, such as convection and diffusion, are not considered in this model, which is acceptable when front propagation is comparatively fast [26]. The appropriate evolution equation for the order parameter is then The decay of radiation intensity as it passes through the mixture is modelled using a generalised Beer-Lambert law given by whereμ is the local optical attenuation coefficient. The dependence of this parameter on φ reflects the fact that the neat liquid phase (φ = 0) and the fully converted phase (φ = 1) may have different attenuation coefficients given by μ 0 and μ ∞ , respectively. To capture this in the model, the local attenuation coefficient is written as the phase average of μ 0 and μ ∞ ; in particular, The model is closed by supplementing it with appropriate boundary and initial conditions. The initial mixture is assumed to be unpolymerized. In addition, the intensity of radiation at the illuminated surface is fixed at I 0 . Therefore, the following conditions are imposed: It is enlightening to first consider the case of photoinvariant photopolymerization, whereby the optical attenuation coefficients of the liquid-and solid-rich phases are equal, that is, μ 0 = μ ∞ ≡μ. In this case, the system of equations in Eq. (1) is readily solved to obtain [16] φ(z,t) The solution for the order parameter can be written in the form of a traveling wave; in particular,  whereẑ is a nondimensional coordinate given bŷ In this case,ẑ measures the distance, normalized byμ −1 , from the inflection point of φ, denoted as z i . Hats are used to distinguish variables that are viewed in a frame of reference that travels with the wave. The profile of this traveling wave is shown in Fig. 2; regions in the "upbeam" and "downbeam" directions,ẑ 0 andẑ 0, correspond to solid-rich and liquid-rich phases, respectively. The wavefront nearẑ = 0 marks the location of a diffuse interface, or interfacial layer, that forms between these two phases. As discussed in Hennessy et al. [26], the nondimensional width of the interfacial layer is given approximately by the inverse of the derivative ofφ at its inflection point:ŵ A geometrical interpretation of this expression is shown in Fig. 2. We find thatŵ = e 1 , which, in dimensional terms, corresponds to w = e 1μ−1 ; thus, under these modeling assumptions, the width is controlled solely by the size of the optical attenuation coefficient. In practice, it is convenient to consider a sharp solidliquid interface that has zero width, defined according to the location where φ reaches a critical value φ c [4,15,18]. This corresponds to the experimentally relevant cutoff of the conversion profile upon development, i.e., selective dissolution of the insufficiently crosslinked liquid-rich phase. We denote this interface by z f = z f (t) and define it implicitly via the expression φ c = φ[z f (t),t]. Using the solution given in Eq. (2a), we find that where In order for (6) to be physically relevant, it must be positive. This requirement is only met if t is greater than the induction time τ ind , i.e., if t > τ ind . This condition reflects the fact that a solid-liquid interface does not instantaneously form when the system is irradiated. Changes in the value of φ c , arising from alternative definitions of the sharp interface based on physical arguments, only modify z f by an additive constant. The fundamental behavior of z f is not influenced by changes in φ c ; therefore, it does not depend critically on the precise definition of the sharp solid-liquid interface. Experimentally, φ c is readily inferred by measuring the position of the planar solidification front z f [4,15] or by resolving the chemical conversion profile φ at a given time t [18]. The model presented in Eq. (1) has been analyzed in the case of μ 0 = μ ∞ by several authors [4,15,16,26]. In general, the solutions of the full model are qualitatively similar to those presented here for the case of photoinvariant photopolymerization.
The remainder of this paper will focus on extending the model in Eq. (1) to include thermal effects in order to examine how these influence the dynamics of FPP, particularly the width of the interfacial layer, w, and the motion of the sharp solidliquid interface, z f .

III. THE ROLE OF THERMAL EFFECTS IN FPP
Extending the minimal model in Eq. (1) to include thermal effects is relatively straightforward. It involves the introduction of an equation governing the temperature field T = T (z,t) and accounting for the temperature dependence of the polymerization rate constant, K. The equation for the temperature is based on conservation of thermal energy and considers three main mechanisms: (i) transfer of heat via thermal diffusion and heat generation due to the (ii) exothermic polymerization reaction and (iii) absorption of incoming radiation. The extended system of equations is given by where ρc p andk are the volumetric heat capacity and thermal conductivity of the mixture, respectively. These material properties could be written as phase averages similar to the local absorption coefficientμ in Eq. (1c); however, for simplicity, they are taken to be constants. The parameter h = H /V m denotes the enthalpy of polymerization per molar volume of monomer and has units of J/m 3 . The temperature-dependent polymerization rate constant K that appears in Eq. (8a) is assumed to have an Arrhenius form, where E a is an activation energy, R is the ideal gas constant, and K 0 is the value of K at a reference temperature of T 0 . The differential equation for the temperature (8b) requires two boundary conditions: one at the illuminated surface and one at the bottom of the mixture. We will assume the system is nonadiabatic so the generated heat can flow out of both boundaries and into the surrounding environment. At the illuminated surface, z = 0, we suppose that the temperature satisfies Newton's law of cooling: where h is a heat transfer coefficient and T 0 is the temperature of the surface. The temperature at the bottom of the mixture is assumed to be fixed at T 0 , which can be controlled, e.g., by placing the system on a hot plate, as regularly done in practice instead of placing the entire lithographic system within an oven. Thus, the second boundary condition can be written as where H is the height of the mixture; see Fig. 1. The initial temperature distribution is given by the steady-state profile in the absence of irradiation: In addition to these equations, the previous boundary and initial conditions for the conversion fraction, φ, and the intensity, I , still apply: We also assume that H is sufficiently large so the bottom of the mixture remains pure in monomer; thus, we suppose that φ(H,t) = 0 for all time. This condition is pertinent for FPP processes, since the pattern height is defined by the front position rather than the thickness of the, e.g., spun cast, photoresist layer. This model is expected to remain valid unless the heat generated by polymerization or radiation absorption is sufficiently large to either (i) activate mass transport mechanisms, such as buoyancy-driven convection and Soret diffusion, or (ii) initiate thermal polymerization. We note here that convection can occur even in the absence of a temperature gradient if the polymer-rich and monomer-rich phases have sufficiently different densities. Neglecting thermal and/or compositional buoyancy-driven convection is an important simplification, as this mode of mass transport can lead to the onset of a Taylor-like instability which, in turn, can drive the breakup of the planar polymerization front [42][43][44].
From this point forward, we will also assume that the local absorption coefficientμ is constant in time and space; thus, we set μ 0 = μ ∞ ≡μ and consider only photoinvariant photopolymerization, which facilitates an analytical study of this model. As shown by Vitale et al. [18], a photoinvariant FPP model can accurately capture experimental data, at least for sufficiently short times. A detailed discussion of the case when μ 0 = μ ∞ is given in Appendix A, where the traveling-wave conversion profile and position of the sharp interface are computed for different ratios of μ 0 and μ ∞ . We find that the temperature influences the dynamics of FPP in similar ways regardless of the relative size of the attenuation coefficients. Therefore, the proceeding discussion, which focuses on photoinvariant photopolymerization, is representative of general photopolymerization.

A. Scaling
The governing equations are now recast into dimensionless form by choosing suitable scales for the variables in the model. This has the advantage of reducing the number of free parameters in the system and it aids in identifying the key parameter and time regimes. We scale the time variable t with (K 0 I 0 ) −1 , which represents the time required for the conversion profile to first have an inflection point in the isothermal case; see (4). In addition, distances are written roughly in terms of the width of the interfacial layer by scaling the space variable z with the inverse of the mean absorption coefficientμ. The temperature is written as the deviation from its initial value of T 0 and is scaled with T , which represents the characteristic increase of the mixture temperature and will be specified below. Finally, we let K 0 and I 0 be characteristic values of the rate constant and intensity, respectively. Using these scales, where primes denote dimensionless quantities. Moreover, we let H =μH be the dimensionless height of the mixture written in terms of the reference interfacial width and assume that H 1. Using the fact that the dimensionless intensity, or the transmission, is given by I (z ,t ) = e −z , the nondimensional model can be written, after dropping the primes, as where α = L diff /μ −1 compares the thermal diffusion length, defined by to the reference interfacial width. The parameter α can also be interpreted as the relative rate of thermal diffusion to polymerization, as faster thermal diffusion would imply a larger value of L diff and hence a larger α. The characteristic temperature increases due to the exothermic polymerization reaction and the absorption of radiation by the mixture, T rxn and T abs , are given by respectively. The nondimensional polymerization rate constant can be written as where are nondimensional numbers. We anticipate [18,45] that the rise in temperature, T , will be small compared to T 0 , which is measured in absolute units (Kelvin). Thus, we have that ε 1, in which case the rate constant in Eq. (16) can be approximated by where ζ = εγ is the Zeldovich number [30]. Unless otherwise stated, we will replace (16) with (18) in the analysis below. The boundary and initial conditions for the nondimensional model are where β is a Biot number given by In the limit β → 0, the illuminated surface becomes a perfect thermal insulator, whereas if β → ∞, then it becomes perfectly conducting. As mentioned earlier, we also assume that φ(H,t) = 0.

B. Analysis
The nonlinear structure of the governing equations prevents an analytical solution from being easily obtained. Thus, we decompose the problem into several simpler problems by making some assumptions about the size of the various nondimensional numbers and hence the physics that drive the evolution of the system.

Recovery of the isothermal model
We first consider one situation where the thermal model reduces to the isothermal model. This occurs when the Zeldovich number ζ is small. In this case, the polymerization constant can be written as K(T ) = 1 + O(ζ ); therefore, to leading order, the reaction rate does not depend on the temperature. Thus, the polymerization reaction evolves independently from the temperature field.

Reaction-dominated heating
We now suppose that the main source of heat generation is due to the exothermic photopolymerization reaction. We still have the freedom to define T so in this regime we choose T = T rxn . To isolate the effects of this heating mechanism, the temperature rise due to radiation absorption is assumed to be small compared to the rise due to the reaction; thus, it is assumed that T abs / T rxn 1. This condition is reasonable if the experiment happens sufficiently rapidly, i.e., if t h/(I 0μ ). The governing system of equations in this regime is where K is given by (18) and with the boundary and initial conditions defined in Eq. (19). a. Limit of slow thermal diffusion. We first consider the limit of slow thermal diffusion relative to the rate of  polymerization. This regime is described by α 1. Our analysis indicates that a thermal boundary layer of width O(α) rapidly develops near the illuminated surface, influencing the local polymerization dynamics. However, we find that this boundary layer does not affect the long-term frontal dynamics so we do not consider it any further. Away from the illuminated surface, i.e., for z α, the diffusive term can be safely neglected from (21b) and the two equations in Eq. (21) can be subtracted to obtain ∂(φ − T )/∂t = 0. This equation can be integrated in time and, by using the fact that T = φ = 0 when z = H , we obtain T = φ for z α. Therefore, the temperature rise follows the polymerization reaction: wherever φ increases, the temperature T does so as well. This solution for the temperature can be inserted into (21a) to obtain The right-hand side of this expression suggests that the reaction can be described by an effective rate constant, This function decreases monotonically if ζ < 1; however, when ζ > 1 it has a global maximum at φ m = 1 − ζ −1 . The solution to (22) can be written as φ(z,t) =φ(ẑ), whereφ is determined from an implicit equation involving exponential integrals [46], andẑ is a traveling-wave coordinate defined bŷ Mathematically, this expression forẑ is equivalent to (4) once the nondimensionalization is taken into consideration; however, the (dimensionless) inflection point of the thermal conversion profile no longer occurs atẑ = 0 as in Eq. (4). The conversion profiles associated with (24) are shown in Fig. 3(a), where they are plotted in such a way that the inflection point ofφ, denoted byẑ i , occurs at the origin. As ζ increases, corresponding to a polymerization rate constant that depends more strongly on the temperature, the conversion profiles sharpen in the sense that the interfacial layer becomes narrower. Recalling that T = φ in this regime, these curves also represent the temperature distribution throughout the mixture. Furthermore, this equality implies that the temperature profile settles into a traveling wave that propagates with the polymerization front. The sharpening of the interface is due to a retention of thermal energy behind the wavefront and the Arrhenius kinetics of the reaction. The retention of heat is caused by the slow rate of thermal diffusion relative to the rate of polymerization, thus preventing the generated heat from spreading ahead of the moving front. The Arrhenius kinetics lead to the conversion of monomer being more rapid behind the front than ahead of it, therefore, giving rise to a sharper, albeit asymmetric, conversion profile. The dependence of the (nondimensional) width of the interfacial layer on the Zeldovich number ζ can be examined by making use of a definition similar to (5), namelŷ The dimensional width of the interfacial layer is w =μ −1ŵ . The inflection point ofφ, as well as the value of the conversion fraction at this point,φ(z i ) ≡φ i , can be determined simultaneously by solving the nonlinear system of equations In general, the solution must be obtained numerically. However, we find that bothẑ i andφ i vary roughly linearly with ζ , suggesting that approximate expressions for these quantities can be obtained by analyzing their behavior for ζ 1 and then extrapolating. Using perturbation methods, the solution of this system of equations for small values of ζ is found to bê By inserting the solution forẑ i into (26), we obtain This approximation breaks down when ζ = (3e −1 − 1) −1 , reflecting the fact that it is technically only valid for ζ 1. Figure 3(b) plots widths that are computed numerically (solid lines) and using the approximation in Eq. (29) (symbols). The two are in excellent agreement, and they both illustrate how moderate increases in ζ lead to large decreases in the width of the interfacial layer.
The (nondimensional) position of the sharp interface defined by φ c and the corresponding induction time can be computed directly from (24) and are found to be The position of the interface retains its logarithmic dependence on time in this thermal regime. Moreover, at each time t, the speed of the interface is equivalent to that found in isothermal systems [see (6)]. However, for a given time, the positions of the interfaces in thermal and isothermal systems will not be the same due to differences in the induction time τ ind . In particular, the accelerated rate of polymerization due to heat release can significantly reduce the induction time, as shown in Fig. 3(c), where τ ind is plotted as a function of ζ for different values of φ c . The monotonic decrease of these curves with ζ reflects the fact that increasing ζ will lead to a greater increase in the reaction rate for a given temperature change, thus enabling the φ = φ c threshold to be crossed in a shorter amount of time. The reduction in the induction time primarily occurs when φ c > 1/2. For φ c 1, the induction time is not strongly altered by thermal effects, even when the coupling between the reaction and the temperature is strong, i.e., when ζ > 1. This is because the interface is formed sufficiently quickly when φ c 1 that strong thermal coupling does not have time to accelerate this process. For ζ 1, we find that τ ind ∼ log[1/(1 − φ c )], thus recovering (7) after the nondimensionalization is taken into consideration.
b. Limit of fast thermal diffusion. We now briefly discuss the case when diffusion of heat is fast compared to the rate of photopolymerization, corresponding to α 1. In this case, the generated heat is rapidly spread across the domain, resulting in small increases in the temperature. It is, therefore, useful to introduce a rescaled temperature θ defined by T (z,t) = α −2 θ (z,t) so T is small when θ = O(1). The governing equations then become subject to We will assume, reasonably, that ζ α 2 . Upon neglecting small terms of O(α −2 ) and higher, we find that the polymerization rate constant no longer depends on the temperature and the leading-order conversion fraction is given by the nondimensional version of the isothermal solution (2). Therefore, the sharpness of the front as well as its propagation remain unchanged in this limit. The temperature evolves quasistatically so at each moment in time it takes on its instantaneous steady-state profile given by where E 1 is an exponential integral that is related to Ei. This solution does not satisfy the initial condition for the temperature, namely θ (z,0) = 0, as it is technically only valid for times which satisfy t α −2 . In fact, there is an initial time regime given by t = O(α −2 ) where the temperature undergoes a rapid transient evolution into the profile given by (32). The solution in this regime can be resolved by first definingt = α 2 t and then letting φ(z,t) = (z,t) and θ (z,t) = (z,t) in the system (31). Upon neglecting terms of O(α −2 ) and higher, we find that (z,t) ≡ 0 and the new temperature variable solves with the following conditions: A series solution to this equation can be obtained; however, it is not required so we do not present it here. It can be shown that (z,t → ∞) = θ (z,t → 0); thus, the solution for the temperature undergoes a smooth transition from the first time regime into the next. c. Effects of moderate thermal diffusion. The dynamics that occur when the time scale of thermal diffusion is commensurate with that of the photopolymerization reaction, i.e., when α = O(1), can be studied by solving (21) numerically. Figure 4 presents the results from a computation when α = 1, β = 0, and ζ = 2, and it shows the evolution of the conversion profile [ Fig. 4(a)], the temperature [ Fig. 4(b)], the width of the interfacial layer [ Fig. 4(c)], and the position of the sharp interface [ Fig. 4(d)]. The heat that is released from the polymerization reaction raises the temperature of the mixture near the illuminated surface. Despite the nonuniformity of the temperature profile, photopolymerization still appears to be frontal, although the conversion profiles are no longer time invariant when viewed from a frame of reference that travels with the interface. As Fig. 4(c) shows, the width of the interfacial layer evolves slowly in time. In particular, the initial width is reduced in comparison to the isothermal case as a consequence of the buildup of temperature behind the wavefront. However, as the system cools due to diffusion, this layer broadens until it reaches the isothermal width. The widths are not shown for t < 1, as for these times the conversion profiles have yet to fully develop an interfacial layer.
The evolution of the sharp interface is shown in Fig. 4 the same. This is due to the interface being rapidly created and propagating into the bulk at a speed that keeps it ahead of the region where the temperature has increased. In the long term, however, this heated thermal region overtakes the interface, with the former advancing like t 1/2 , due to diffusion being the driving factor, and the latter like log t. Once the interface is overtaken, it experiences an acceleration due to increased rates of polymerization and generation of solid-rich phase. As the system cools, the interface decelerates and begins to travel at the same speed as in the isothermal case. What is particularly remarkable is that the gain in position from thermally induced acceleration is lost during the deceleration phase and, in the long term, the interface evolves as if the system had been isothermal the entire time. This return to isothermal behavior is most easily explained if Fig. 4(d) is interpreted as showing the time t at which an interface forms at a distance z f from the illuminated upper surface. Keeping this in mind, we find that the temperature rise decreases with z f ; therefore, the polymerization kinetics are not subject to a strong thermal acceleration far from the upper surface. Hence, the interfaces should form at roughly the same time as in the isothermal case at these positions. Similar behavior is seen for larger values of φ c , although the acceleration of the interfacial position happens immediately due to thermal effects influencing the polymerization reaction before φ c is reached.
From numerical experimentation we have found that the qualitative aspects of these results remain unchanged under parameter variation. The most significant change occurs when β increases. When this is the case, more heat is able to flow out of the system through the illuminated surface. This reduces the temperature in the mixture and hence also the interfacial sharpening and acceleration.

Absorption-dominated heating
We now consider the opposite physical situation where the dominant heating mechanism is radiation absorption rather than exothermic reactions. Thus, we set T = T abs and suppose that T rxn / T abs 1. The governing system of equations in this case can be written as with K as in Eq. (18) and boundary and initial conditions given by (19). The problem for the temperature decouples from that for the conversion fraction and it can be solved analytically using a Fourier series. By solving (34a), we find that the conversion fraction can be written in terms of the temperature: Although a closed-form expression for the temperature is known, it is sufficiently complicated to prevent the integral in Eq. (35) from being evaluated analytically. a. Limit of slow thermal diffusion. An interesting solution can, however, be found in the limit of slow thermal diffusion, i.e., when α 1. We find that thermal boundary layers form at both the top and bottom of the mixture, the former being analogous to that seen in the case of reaction-dominated heating and the latter only being noticeable if the mixture is not very deep. These thin layers do not affect the structure of the travelling waves; hence, we do not consider them further. In the interior of the mixture, the diffusive term can be neglected from (34b) and the simplified equation can be solved to find T (z,t) t exp(−z). Thus, the temperature of the mixture increases linearly with time from the constant absorption of radiation but decays exponentially from the surface due to the spatial attenuation of light. In fact, the temperature can be written as T (z,t) = tI (z), showing that the temperature is simply a scalar multiple of the intensity profile. By inserting the solution for the temperature into (35), we find that whereẑ = z − log t. Strictly speaking, this solution is only valid for nondimensional times which satisfy t ε −1 , where ε is given in Eq. (17). For larger times, the temperature is sufficiently large that the full form of K given by (16) must be used. The inflection point of (36) satisfies the nonlinear equation and once it is determined, the width of the interfacial layer can be evaluated using (26). The location of the sharp interface and the corresponding induction time can be determined directly from (36) and are found to be z f (t) = log(t/τ ind ), (38a) The dynamics of the system in this regime are qualitatively very similar to those in the case of radiation-dominated heating with slow thermal diffusion. Increasing ζ leads to conversion profiles with thinner interfacial layers, as shown in Fig. 5(a). The interfacial sharpening mechanism parallels that seen previously as well. In particular, the temperature can be written as T (ẑ) = e −ẑ , which implies that it remains constant in time when viewed from a frame of reference that travels with the polymerization front. The slow rate of thermal diffusion maintains high temperatures behind the front and low temperatures ahead of it. The polymerization reaction is, therefore, accelerated behind the interface, whereas it proceeds at essentially the same rate as in the isothermal case ahead of it.
The width of the interfacial layer as a function of ζ is shown in Fig. 5(b). The same qualitative behavior is observed here as in the case of reaction-dominated heating, although the dependence ofŵ on ζ is now much weaker. In addition, the inflection pointẑ i varies nonlinearly with ζ , making it difficult to obtain an approximate expression for this quantity that is uniformly valid over a wide range of values. However, it is possible to obtain expressions that are valid for small and large ζ . For ζ 1 we find that whereas for ζ 1, These approximations are compared to numerical solutions in Fig. 5(b) and the agreement is found to be excellent. Examining (38) shows that the position of the sharp interface again has a logarithmic dependence on time. Moreover, the variation of the induction time with ζ and φ c , shown in Fig. 5(c), follows the same trends as in the case of radiation-dominated heating.
b. Limit of fast thermal diffusion. The analysis of the model in the case of fast thermal diffusion, α 1, follows the same approach as before. We introduce a new temperature variable θ , defined according to T (z,t) = α −2 θ (z,t), and suppose that t α −2 and ζ α 2 . After neglecting terms of O(α −2 ) and higher, we find that the leading-order solution for the conversion profile is the same as in the isothermal case and given by (2). The temperature is in a steady state which can be written as The solution in the t = O(α −2 ) time regime satisfies the same problem as in Eq. (33). This implies that in the limit of fast thermal diffusion, the only way to distinguish a system that is heated primarily by the exothermic polymerization reaction rather than the absorption of radiation is through the longterm behavior of the temperature field. In the former case, the temperature rise at a given point will decay to zero, whereas in the latter it will reach a constant nonzero value.
c. Effects of moderate thermal diffusion. As before, the case when α = O(1) can be studied via numerical simulations of the governing equations (34). Figure 6 shows the results of a simulation using α = 1, ζ = 1, β = 1, H = 30, and φ c = 0.05. The transient dynamics are very similar to those seen when the heating is due to the polymerization reaction and will not, therefore, be discussed further. The long-term dynamics differ considerably, however, due to the nontrivial steady-state temperature profile, shown in Fig. 6(b). As a result, the conversion profiles are permanently sharpened [ Fig. 6(c)] and for a given time, solid-liquid interface has always propagated further into the bulk than in the isothermal case [ Fig. 6(d)]. The solid-liquid interface experiences a marked acceleration that begins at t 10, corresponding to the time when it is overtaken by the heated region of the mixture. For larger times that roughly satisfy t > 10 3 , the interface decelerates slightly as it propagates into a colder bath, reducing the gain it had over the isothermal interface.
To further explore the behavior of the system at large times, we compute asymptotic expressions for the conversion profile, the width of the interfacial layer, and the position of the interface by considering the evolution after the temperature has approximately settled into its steady-state configuration. We define a time t 1 1 such that T (z,t) α −2 θ ∞ (z) for t > t 1 , where θ ∞ is given by (40). In addition, we let φ(z,t 1 ) = φ 1 (z). In practice, φ 1 can be determined by solving the transient problem starting from t = 0; however, a precise functional form is not required here. If we consider only the linear contribution to the temperature profile, which turns out to be the only relevant contribution for large times, and neglect exponentially small terms proportional to e −H , then it can be shown that the order parameter satisfies with φ(z, for t t 1 1. The width of the interfacial layer and the position of the sharp interface are given by which hold for t t 1 1. Using the expression in Eq. (42c), it can be shown that the sharp interface reaches the bottom of the mixture at the same time as in the isothermal case. This happens when t = e H log[1/(1 − φ c )], independent of ζ . Therefore, the long-term kinetics of wave propagation are ultimately enslaved to the "cold" boundary at z = H , which serves to constantly diminish any gains that are obtained through thermal effects, regardless of the coupling strength between temperature and the reaction rate.
The asymptotic approximations for the conversion fraction (42a), width of the interfacial layer (42b), and the position of the sharp interface (42c) are shown as symbols in Figs. 6(a), 6(c), and 6(d), respectively. The approximations clearly capture the long-term behavior of the numerical solution.

IV. DISCUSSION AND CONCLUSION
In this paper we have systematically investigated how the generation of heat within a photopolymerizing bath can influence the onset and dynamics of FPP. Heat generation was assumed to be due to exothermic reactions and the absorption of radiation. Particular attention was paid to determining the shape of traveling-wave conversion profiles and the motion of a sharp solid-liquid interface defined by the location where the conversion fraction reaches a critical value. Both of these aspects are directly relevant to the practical use of FPP. For instance, the position of the sharp interface determines the patterned size of the growing network solid. The ability to control the width of the interfacial layer can facilitate the fabrication of gradient polymer solids with material properties that vary in a specific manner.
Heat generation during photopolymerization creates a thermal front that propagates into the bath with the polymerization front. Depending on the relative rates of thermal diffusion and photopolymerization, the thermal front can advance as fast, or faster, than the polymerization front. We find, in particular, that three key regimes emerge. When thermal diffusion is fast compared to photopolymerization, the thermal front propagates much more rapidly than the polymerization front [ Fig. 7(a)]. The fast distribution of heat across the large domain results in small temperature increases that do not significantly affect the dynamics of FPP. If diffusion is comparatively slow, the thermal and polymerization fronts travel together at the same rate [ Fig. 7(b)]. The accumulation of thermal energy behind the polymerization front in combination with Arrhenius reaction kinetics leads to a sharpening of the conversion profile. Moreover, as the generated heat cannot spread ahead of the polymerization front into monomer-rich fluid, the solid-liquid interface evolves as if the system were isothermal. Finally, when the time scales of diffusion and photopolymerization are similar and φ c is small [ Fig. 7(c)], the initial propagation of the thermal front will be as fast as the polymerization front. However, the thermal front, being driven by diffusion in the long term, eventually overtakes the polymerization front, at which point the evolution of the latter becomes nonlogarithmic. Remarkably, the logarithmic behavior of the polymerization front is recovered for large times due to the system cooling to its initial temperature or the temperature settling into a linear profile. For larger values of φ c , the thermal front is always ahead of the polymerization front due to increased induction times. We found that FPP can occur in all three of these regimes, although the conversion profiles and propagation rates can differ slightly from their isothermal counterparts.
An interesting prediction of our model is that a buildup of heat behind the polymerization front can influence the shape of the conversion profile. This was demonstrated particularly clearly when thermal diffusion and polymerization occur on similar time scales and the generation of heat is primarily due to radiation absorption; in this case, the presence of a strong and permanent temperature gradient led to long-term changes in the width of the interfacial layer as well as the motion of the sharp solid-liquid interface. Both of these quantities were found to depend on a single nondimensional number , which can be written in terms of dimensional quantities as The direct relationship between and dT /dz suggests that the dynamics of FPP can be tailored by imposing a temperature gradient along the growth direction of the solid, therefore, enabling the coupling between temperature and photopolymerization to be taken advantage of in a precise manner. We will explore this novel approach to controlling FPP in an upcoming publication.
A key assumption underlying the thermal model is that the temperature increases are sufficiently small that mass transport, e.g., convection and diffusion, can be neglected. In the limit of fast thermal diffusion, the temperature increases are always small so this assumption is likely to hold true.
When thermal diffusion is slow, however, there are long-lived temperature increases which could be sufficiently large to activate mass-transport mechanisms. If this were the case, then the sharpening of the diffuse interface due to thermal effects would compete with broadening of it due, for example, to mass diffusion [26]. A beneficial future study would, therefore, be to further ascertain the relationship between the rate of thermal diffusion, the onset of mass transport, and the dynamics of FPP. Further, since many light-driven polymerization processes can also be thermally initiated if a thermal-sensitive initiator is added to the reactive mixture, we expect the coupling of light-and heat-induced polymerization to have important implications for manufacturing.
The thermal model presented here implicitly assumes that diffusion of heat is restricted to one spatial dimension, which strictly applies when the lateral dimensions of the monomer bath are much greater than its height. In the case of lateral confinement, e.g., when the photomask aperture is of same order as the height of the bath, there can be a significant transfer of heat to the cooler liquid-rich phase through the sides of the solid that are orthogonal to the illuminated surface. The resulting lateral temperature gradient, with the center of the solid being warmer than its sides, in combination with Arrhenius reaction kinetics, could drive the formation of a curved solidification front. However, such a scenario is only likely to occur in the case when the dimensionless parameter α is order one in size. For large values of α, diffusion of heat will still be sufficiently fast to keep the mixture approximately uniform in temperature; when α is small, diffusion of heat through the sidewalls of the solid will be too slow to reduce the local temperature and hence impact the process of photopolymerization. Thus, in both the large-α and small-α limits, the dynamics of FPP are essentially one-dimensional, regardless of the lateral size of the solid.
Applying the thermal model in practice requires a careful consideration of the boundary conditions at the bottom of the mixture. An underlying assumption of our analysis is that the bath is sufficiently deep for the temperature at the bottom to remain constant and equal to its initial value. For shallower baths, such as those where the depth is on the order of the interfacial width, i.e., H = O(μ −1 ), there can be a significant accumulation of heat near the bottom boundary due to generation or diffusive transport. This, in turn, can lead to a local increase in temperature if this boundary is at least partially insulating. Thus, when using the model to describe shallower mixtures, we recommend replacing the boundary condition (10) with an analogous version of (9), which would account for such behavior. Indeed, we have shown that by modifying the thermal model in this way, it can accurately capture experimental data obtained in several nonisothermal FPP systems [45].
This paper highlights how the additional mathematical tractability afforded by a simple, minimal model of photopolymerization can yield novel insights that are relevant for practical uses of this process. In addition, the ability to present our theoretical results in terms of experimentally measurable quantities will facilitate comparisons between the model and reality, therefore leading to a better quantitative understanding of FPP and how it can be used to fabricate state-of-the-art polymer network solids.

ACKNOWLEDGMENTS
We thank the Engineering and Physical Sciences Research Council (EPSRC) for "Manufacturing with Light" grant No. EP/L022176/1 and Jack. F. Douglas (NIST) for stimulating discussions.

APPENDIX A: OTHER PHOTOPOLYMERIZATION TYPES
Section III focused on the case of photoinvariant photopolymerization, whereby the optical attenuation coefficients of the solid-and liquid-rich phases are equal. We now generalize the analysis to the case when these coefficients are not equal and show that the key results also apply to more general FPP situations.
It is straightforward to show that in the limit of fast thermal diffusion, α 1, the leading-order problem reduces to that of the minimal model of Sec. II for all combinations of μ 0 and μ ∞ . Therefore, the system will exhibit isothermal behavior.
Extending the analysis in the limit of slow thermal diffusion, α 1, is more technical and will be the focus of the remainder of this appendix. In general, the analysis cannot be carried out analytically, although substantial progress can be made in the limits of strong photodarkening (μ 0 μ ∞ ) and strong photobleaching (μ ∞ μ 0 ). We find that the width of the interfacial layer is set by the dominant optical attenuation coefficient and is a decreasing function of the Zeldovich number ζ . In addition, the temperature still follows the conversion and intensity profiles in the cases of reaction-and absorption-driven heating, respectively, and the scaling laws for the position of the interface with time remain the same.

Strong photodarkening and slow thermal diffusion
We first consider photodarkening photopolymerization by settingμ(φ) = μ ∞ φ in the thermal model (8). The equations are then nondimensionalised as in Sec. III A withμ replaced by μ ∞ . In the case of reaction-dominated heating, i.e., T abs / T rxn 1, the governing equations become with φ = T = 0 at t = 0 and I = 1 at z = 0. We will assume, for simplicity, that the domain is unbounded and that φ,T ,I → 0 as z → ∞. We seek a traveling-wave solution to these equations of the form the system reduces to Upon dividing these equations, we obtain subject toÎ (φ = 0) = 0. This equation can be solved to find that By inserting this expression into (A3a), a single differential equation is obtained for the conversion profile. To close the problem, the conditionφ(0) = φ c is imposed. Solutions to this equation are shown in Fig. 8(a).
In the case of absorption-dominated heating, T rxn / T abs 1, the leading-order equations are given by with φ = T = 0 at t = 0, I = 1 at z = 0, and φ,T ,I → 0 as z → ∞. We again seek a traveling-wave solution of the form (A2c) and immediately find thatT =Î . Thus, the problem reduces to Following the same procedure as above, we find that thus leading to withφ(0) = φ c . Solutions of this equation are shown in Fig. 8(b).

Strong photobleaching and slow thermal diffusion
In this case we setμ(φ) = μ 0 (1 − φ) in the thermal model (8) and nondimensionalize using μ 0 rather thanμ. The analysis proceeds in exactly the same manner as in Appendix A 1 so we omit the details here. In the limit of reaction-dominated heating, we find that withφ(0) = φ c . The conversion profiles in this case are shown in Fig. 8(c).
In the limit of absorption-dominated heating, the solution satisfiesT subject toφ(0) = φ c . Figure 8(d) shows the conversion profiles in this case.

Determination of the interface position
for radiation-dominated heating, anḋ for absorption-dominated heating. In both cases, the sharp interface propagates at a constant speed throughout the mixture at a rate that increases with ζ . In addition, the long-term behavior is independent of φ c , indicating that thermal effects will eventually become important, even if φ c 1. This is in stark contrast to what is found for photoinvariant photopolymerization, where the velocity is given byż f = t −1 , which not only decreases with time but is also independent of ζ .
Obtaining the long-term behavior of the interface in the case of strong photodarkening is more involved, as the intensityÎ is singular whenφ = 1 and, therefore, it cannot simply be expanded about small δ. Instead, a precise functional form for δ = δ(ẑ), valid forẑ 0, will be determined and used to construct the leading-order part of (A17), which then can be integrated. In particular, we letφ(ẑ) = 1 − δ(ẑ) in Eq. (A3a) with (A5), or (A9), and then expand about δ 1. For reactiondominated heating, we find that whereas for absorption-dominated heating, where A > 0 is a constant of integration and g is a complicated function of ζ with the property that g(0) = 0. In both cases, the leading-order part of (A17) is given byż f ∼ B(ζ )e −z f , where B is a constant depending on ζ . The long-term motion of the interface is, therefore, logarithmic and given by z f (t) ∼ log t + log B(ζ ), t 1.
Differentiating this expression with respect to time shows the the velocity of the interface is independent of ζ , similarly to the case of photoinvariant photopolymerization. However, the time at which the transition to logarithmic behavior occurs does depend on ζ . The asymptotic expressions for the position of the sharp interface are validated using numerical simulations of the thermal model in the limit of slow thermal diffusion. Figure 9 shows the evolution of this interface in the cases of strong photodarkening and strong photobleaching with radiation-and absorption-dominated heating. In all cases, we have set ζ = 2 and φ c = 0.05. The asymptotic solutions (dashed lines) are in excellent agreement with the numerics (solid lines). Moreover, Figs. 9(a) and 9(b) clearly show the transition from algebraic behavior, given by (A15b), to logarithmic behavior predicted for photodarkening photopolymerization.

APPENDIX B: THE UNIVERSAL LARGE-TIME LIMIT OF THE MINIMAL MODEL
A remarkable property of the minimal model in the case of photoinvariant photopolymerization is that the conversion fraction approaches the same travelling-wave profile regardless of the initial condition. This implies that the initial evolution of the system is eventually forgotten about. To see this in detail, consider the initial value problem given by subject to φ(z,t 1 ) = φ 1 (z), where A and B are arbitrary constants. It is reasonable to assume that φ 1 satisfies φ 1 (z → ∞) → 0, so no polymerization has taken place in the far field. The solution of this problem is given by We rewrite this solution in terms of the coordinateẑ = z−s(t), where to obtain At the moment, s has no physical meaning, as it corresponds to neither the inflection point of the conversion profile nor the point where a critical value of φ c is obtained. By taking the limit t t 1 , we find that Therefore, in the long-time limit, s converges to the inflection point and the conversion profile becomes independent of the initial condition.