Modeling of Photochemical Reactions in a Focused Laser Beam

Fluorescent materials play a prominent role in the qualitative and quantitative measurement of scientific phenomena of importance in biotechnology and biomedical applications. Photodegradation of fluorophores is a process that determines the accuracy and sensitivity of such measurements. This is the motivation for developing methods for accurately measuring fluorophore photodegradation rates. Recently, illumination consisting of short pulses has been used to examine the decay of photochemical reaction products. However, the time resolved measurements are difficult to interpret since the photodegradation process usually involves multiple time scales. The frequency domain measurement technique discussed here looks at the frequency response of a fluorescent sample to a frequency modulated illuminating light. The photodegradation rate is obtained by interpreting the frequency domain measurements in terms of traditional impedance concepts. In the measurements described in this paper, a focused laser beam is used to illuminate a sample of slowly flowing fluorescent solution. The laser beam is assumed to have a Gaussian power distribution hence illumination is spatially non-uniform in the region of interest. The photochemical reaction rates depend on power, so they will also vary with the position in the beam. However in the case of photodegradation of fluorophores, the measurement of the resulting decrease in fluorescence is given in terms of the radiation emitted from the entire illuminated region. In this work we present a mathematical description of the time evolution of the fluorescence response integrated over a non-uniformly illuminated domain. As a result of our analysis, an experimentally accessible and tractable mathematical model Eq. (19) and Eq. (30) is obtained from a more fundamental description given by Eq. (4) and Eq. (5). The model is used to create a functional form for fitting experimental measurements from a lock-in amplifier.

Fluorescent materials play a prominent role in the qualitative and quantitative measurement of scientific phenomena of importance in biotechnology and biomedical applications. Photodegradation of fluorophores is a process that determines the accuracy and sensitivity of such measurements.This is the motivation for developing methods for accurately measuring fluorophore photodegradation rates. Recently, illumination consisting of short pulses has been used to examine the decay of photochemical reaction products. However, the time resolved measurements are difficult to interpret since the photodegradation process usually involves multiple time scales. The frequency domain measurement technique discussed here looks at the frequency response of a fluorescent sample to a frequency modulated illuminating light. The photodegradation rate is obtained by interpreting the frequency domain measurements in terms of traditional impedance concepts. In the measurements described in this paper, a focused laser beam is used to illuminate a sample of slowly flowing fluorescent solution. The laser beam is assumed to have a Gaussian power distribution hence illumination is spatially non-uniform in the region of interest. The photochemical reaction rates depend on power, so they will also vary with the position in the beam. However in the case of photodegradation of fluorophores, the measurement of the resulting decrease in fluorescence is given in terms of the radiation emitted from the entire illuminated region. In this work we present a mathematical description of the time evolution of the fluorescence response integrated over a non-uniformly illuminated domain. As a result of our analysis, an experimentally accessible and tractable mathematical model Eq. (19) and Eq. (30) is obtained from a more fundamental description given by Eq. (4) and Eq. (5). The model is used to create a functional form for fitting experimental measurements from a lock-in amplifier.
below, a focused laser beam is used to illuminate the sample. The power distribution of a Gaussian laser beam is strongly dependent on spatial location; therefore, the rates of photochemical reactions, which depend on power, will also depend on location in the beam. Since absorption of radiation leads to local heating and concentration gradients, convective and diffusive mass transport will be present in the vicinity of the focused beam. A slow flow is imposed on the fluorophore solution in order to dominate and thus minimize the effect due to uncontrollable convective and diffusive mass transport. In the case of photodegradation of fluorophores, measurement of the resulting decrease in fluorescence usually measures the radiation emitted from the entire illuminated region. Therefore an analysis method has to be developed which describes the time evolution of the integrated fluorescence response from a nonuniformly illuminated region. In previous work [6] a frequency domain measurement technique was interpreted with a model which assumed a uniform beam profile. This work extends the model to a laser beam with a Gaussian profile. In it we derive an experimentally accessible mathematical model Eq. (19) and Eq. (30) from first principles Eqs. (4) and (5). A dramatic simplification of these equations is a consequence of the strongly disparate time scales in the kinetics imposed by the slow rate of degradation relative to the relaxation and absorption rates The first section contains a motivation and presentation of the equations that are taken as a fundamental description of the experimental setting (also see Fig. 1). The equations are based on a two-state model with excitation and relaxation. In Sec. 2, a perturbation scheme employing the disparate time scales is used to reduce the system of two partial differential equations, describing the kinetic process to a single partial differential equation Eq. (14) for the total fluorophore population (to leading order in the perturbation parameter). The biggest reduction in mathematical complexity comes from writing the total measured fluorescence in terms of the total fluorophore population averaged over the Gaussian light beam distribution. This formula is introduced in Sec. 3, Eq. (9) and there follows the derivation of a single ordinary equation for the time evolution of the spatial average Eq. (30) that completes the model. The mathematical arguments that justify the manipulations performed in the main body of the paper can be found in the appendices. Section 4 discusses the time independent case and a functional form for fitting experimental measurements from a lock-in amplifier. The form, obtained from fitting values predicted from the model, shows good agreement with actual data. In a forthcoming paper, a harmonic approximation of the solution of the model equations will be developed. It will be possible then to express the fit parameters in terms of the physical parameters in the experiment.

Kinetic Equations
As a starting point, consider a model for reactions at a fixed point in space. This simplified model will be used to motivate the details of the analysis. Assuming no mass transfer from surrounding fluid, Fig. 1 shows a two state model where absorption occurs from the ground state with rate k a , and the excited state population of size N 1 , relaxes back to the ground state with a rate k r . Here N 0 is the population of the ground state molecules. Intact fluorophores change into non-fluorescent species with a rate k d . The kinetics are summarized by a set of equations (1) To interpret k d , Eq. (1) has to be supplemented by a model of the photochemical processes that are responsible for the conversion of the fluorophores into nonfluorescent species [6]. In the present manuscript, there will be no attempt to interpret the rate constant k d . The discussion below will focus on two problems inherent Fig. 1. A graphical representation of the components included in the kinetic model discussed in the text. The rate constant k d , represents photodegradation that leads to non-fluorescent species. The rate constants k a and k r represent absorption and relaxation to the ground state respectively. N 0 and N 1 represent the populations of the ground and excited states.
in using Eq. (1) to describe real experimental conditions. The first problem is that the laser beam has a spatial intensity distribution leading to spatial variation of the absorption rate. The second problem is that convection and diffusion currents will change the local concentration of molecules in an uncontrolled manner.
The second problem will be reduced significantly by introducing a flow that will change the concentration of fluorophores in a predictable way. The absorption rate depends on the power density of the incident light. The absorption rate will be written as k a = σ a I [7], where σ a (cm 2 ) is the molecular absorption cross section, and I is the incident photon flux (1/s cm 2 ). During the experiment, the laser irradiance P(W/cm 2 ) is measured. The irradiance P can be converted into a photon flux, I, by dividing P by the energy per photon. Explicitly, I = P c P where the conversion n is the index of refraction, h is Planck's constant, and c is the speed of light in vacuum. This yields the result k a = σ a P c P = k c P. The constant k c relates the laser irradiance to the molecular absorption rate. Figure 2 shows the geometric relation of the laser beam, the cuvette holding the sample, the flow, and the detector. The laser beam is incident along the z axis which points out of the plane of the paper. The irradiance of the modulated beam that illuminates the sample will be written as: (2) P 0 (x,y) is the time independent component, and ∆P(x,y) is the amplitude of the modulated component, ω is the and t is time. The beam will be assumed to have a Gaussian profile [8].
where P 0 is the total power (watts) of the laser beam and w is the width of the beam. The function f (x,y) is defined such that its integral over all the x,y plane is normalized to 1. The time independent component and the modulation amplitude have the same spatial depend-ence given by f (x,y). Clearly in this situation, the absorption rate will have a strong spatial dependence. Consequently reaction products may have concentration gradients and uneven heating may cause thermal mass transfer. The kinetic model of photodegradation, given by Eq. (1), will be extended to include a spatially dependent absorption rate and a constant flow of the solution containing the reaction product. It will be assumed that the flow dominates all mass transport (diffusive and convective) into and out of the illuminating region and provides a well defined initial concentration of fluorophore in the illuminated region.
To account for the spatial variation of absorption and the flow of fluid carrying fluorophores past the laser light, Eq. (1) is modified to,  Here k a (x,y,t) = k c P(x,y,t). The velocity of the fluid, denoted by v is assumed to be a constant and in the direction of the x axis in Fig 2. Note in contrast to Eq. (1) the absorption rate is now a function of space and time. To completely specify the solution of these equations, the number of particles flowing into the entrance of the apparatus must be specified and the number of particles at the beginning of the observation period (t = 0) must also be given. Equation (5) states that prior to entering the laser beam there is no change in the populations of the fluorophore states. The distance L is some length along the x axis where the laser beam intensity is effectively zero. Eq. (5) also gives the initial conditions at t = 0 when the laser beam is turned on.

Perturbation Analysis of Kinetic Equations
The faster processes, absorption and relaxation dominate the time variation at least initially. Returning to the model for the fixed position in Eq. (1), this suggests that after an initial transient the quasi-equilibrium condition holds, Substitution of this condition into Eq. (1) results in a single equation for the total fluorophore population N = N 0 + N 1 .
In this section this reduction is extended to the case of non-uniform illumination and is made rigorous. Here the coupled pair of partial differential equations, Eq. (4), will reduce to a single partial differential equation Eq. (14), accurate to first order in a perturbation expansion. Justification of the operations leading to these simplifications can be found in Appendix A.
The absorption and the excited state relaxation rates are usually of the order 10 8 s -1 , while the value of k d is less than 10 3 s -1 . The constants in Eq. (4) can therefore be re-scaled as k a = k a / ε, k r = k r / ε where ε = 10 -8 and k a , k r are of the order of 1. Equation (4) can then be rewritten as, We will treat N 0 (x,y,t,) and N 1 (x,y,t) as functions of ε. Let N ∧ 0 (x,y,t) = N 0 (x,y,t,ε) | ε = 0 and N ∧ 1 (x,y,t) = N 1 (x,y,t,ε) | ε = 0 be the leading terms for N 0 (x,y,t,ε) and N 1 (x,y,t,ε) respectively in a perturbation expansion for ε > 0 and large time t.
The validity of the expansion in ε is discussed further in Appendix A. Substituting Eq. (7) into Eq. (6) and equating coefficients of equal powers of ε on both sides of the equation gives for ε = 0 (8) This implies that as ε → 0 (9) Thus the relation between the populations of ground and excited states holds in this approximate sense for all larger times. (This is not valid for t ≈ 0. When the illumination is turned on, there is an initial transient period where the populations are changing rapidly and therefore where the left hand side of Eq. (6) is not proportionate to ε.) Equating the coefficients of the first power of ε gives: , be the total fluorophore population in the zero order approximation, we can add the two parts of Eq. (10) and obtain (11) The right side of Eq. (11) depends on the zeroth order approximation to the population of excited state. Using Eq. (9) and the definition of the total fluorophore population it is possible to obtain the relation (12) which is correct to zero order in ε. Substituting Eq. (12) into the right hand side of Eq. (11) gives the following equation for N (x,y,t) valid up to first order in ε The boundary condition for N(x,y,t) at x = -L is inherited from those imposed on N 0 and N 1 at The initial condition at t = 0 is N(x,y,0) = N 0 since no reaction takes place prior to turning on the illumination.
The last term in Eq. (14) gives the flow related flux of fluorophores into the illuminating region as a function of position. The flow of fluorophores into the volume illuminated by the laser beam is assumed to be uniform over the entire volume. The photodegradation depends on the intensity of incident light; consequently the concentration of intact fluorophore will not be uniform over the illuminated volume.

Spatial Averaging and the Fluorescent Signal
The measured fluorescence signal per unit distance in the z direction is given by (16) where A represents the characteristics of the optical measurement instrument, k rad is the rate of radiative decay of the excited state population, and N 1 (x,y,t) is the population of the excited state. (The radiative decay contributes to the overall relaxation rate given by k r ). The population of the excited state can be approximated by Eq. (12) leading to (17) The symbol b was defined previously. Equation (17) together with Eq. (14) constitute a description of the measurement valid to first order in ε. Equation (17) suggests that instead of finding the solution N(x,y,t) at every spatial point it may be better to find the variation of the averaged concentration defined by Eq. (18) The average is performed with a weighting factor equal to the spatial distribution of the laser irradiance. The weighting factor provides a measure of the contribution of fluorophores at each point in space to the total fluorescence signal. It will be shown below that Eq. (14) can be converted into an equation for 〈N(t)〉 as defined by Eq. (18), and that the fluorescence signal can be written as where to a very good approximation α is a constant and P(t) = P 0 + P 1 cos(ω t) ≡ P 0 + ∆(t). For actual computations, the constant α is evaluated in terms of the time independent solution of Eq. 14. In this section, an ordinary differential equation for the time evolution of 〈N(t)〉 will be derived (Eq. (30)). The fluorescent signal can therefore be obtained (with small error) by solving Eq. (30) and then using Eq. (19). Together these equations represent an experimentally accessible mathematical model of the fluorescent signal under excitation by a focused laser beam.
To find an equation for 〈N(t)〉, first multiply Eq. (14) by 1 + bP(x,y,t) and obtain (below r represents the coordinates x,y). (20) Multiply Eq. (20) by f (x,y) and integrate over the beam cross section region, R. Each term will be treated in detail. To simplify the notation integrals in the x-y plane will be represented by a single integral in r = (x,y) with dr = dxdy. The left side of Eq. (20) becomes which can be expanded to Inserting the explicit form of the power function we get The second part of the above equation will be rewritten by introducing the quantity, The final result for the left side of Eq. (20) after multiplication and integration is therefore (22) (Note that the definition given in Eq. (21) will be used in deriving Eq. (19)).
The first term on the right side of Eq. (20) after the operations becomes (23) where ais defined in Eq. (21). The second term on the right side of Eq. (20) becomes, after multiplication and integration, To treat the second term of the above equation, we introduce a second function defined by (24) The second term of the right hand side of the transformed Eq. (20) now reads (25) In Appendix B we will show that the functions α -(t) and α -′(t) are, up to a small error, independent of t and can be replaced by constants α and α ′ respectively. The error introduced into equations Eq.
Collecting terms, Eq. (20) can be written after multiplication by f (x,y) and integration, up to small varying functions of t, as where (27) The final step in the derivation of an equation for 〈N(t)〉 consists of finding an approximation to the weighted average spatial derivative in Eq. (27). For small widths w the major contribution to the integral comes from a region in the vicinity of the beam which we define by 0 ≤x,y≤ κ w where we may choose the constant κ , 0 ≤ κ ≤ 3. If we replace the partial derivative in the integrand by a backwards difference quotient with step h =κ w the value of N at the left endpoint of the interval is N 0 (to a good approximation). Thus for w small enough we have, where o(1)is a term that tends to zero as w → 0. To guarantee an error that is smaller than the dominant terms of Eq. (26) we have to require that w is sufficiently small and κ is large. In fact experimental conditions govern the actual choice of w. The following heuristic argument governed the choice of κ.
We observed that in the vicinity of the origin (x = 0, y = 0 ) the spatial derivative of the time independent solution, N(x,y) discussed below, has a functional form which is almost identical to the beam distribution function f (x,y) as shown by the two overlapping curves in Fig. 4. The difference between -f (x,y) and ∂N(x,y)/∂x (both normalized by dividing by the maximum of their respective absolute values) is magnified 50 times and shown in Fig. 4 by the upper curve. The weighted average of the un-normalized beam distribution function π ⋅ width 2 f (x,y) is equal to 0.5. Therefore the weighted average of the spatial derivative will be approximately equal to 0.5 times the value at the origin (28) The derivative at the origin can be approximated in the usual way by the difference of the function N(x,y,t) evaluated at two points spanning the origin and divided by the distance between the two points. However for the purpose of this analysis a more useful approximation uses the difference 〈N (t)〉 -N 0 divided by an effective distance which reproduces the value of the derivative. For the case of the Gaussian beam, the best estimate of the derivative at the origin was found (using the time independent solution) to be (〈N (t)〉 -N 0 )/w*0.637. This yields the following approximate expression for the average spatial derivative given in Eq. (28).

Validation of the Model Equations and Their Derivation:
This section discusses the time independent form of Eq.
(31) Figure 2 shows the geometrical arrangement of the flow channel, the laser beam, and the detector. The laser beam is along the z direction, the detector along the y axis, and the flow takes place along the x axis. The channel dimensions are d y = 0.4 cm along the y axis, d z = 1 cm along the z axis, and 5 cm along the x axis. The velocity of the fluid in the center of the channel is set to 0.025 cm/s. No uncertainties are given since these values will be used to model the response of the system. The width of the Gaussian laser beam is set to w = 50 ⋅ 10 -4 cm. The laser power is P 0 = 0.03 W, the vacuum wavelength is λ = 488 ⋅ 10 -7 cm. The number of photons per Joule, P c , is calculated by P c = λ n/hc = 3.314 ⋅ 10 21 J -1 , where n = 1.35 is the index of refraction of water and hc is the product of Planck constant and the speed of light. The extinction coefficient is set to ε = 82000 L/cm mol, yielding a cross section σ = 1000ε /0.4343 N A = 3.135 ⋅ 10 -16 cm 2 where N A is the Avogadro's number. The average absorption rate is given by k a = σ P c P 0 /√π width = 3.97 ⋅ 10 5 s -1 , and the relaxation rate is given by the inverse of the measured lifetime of the excited state, k r = 1 /τ = 3.33 ⋅ 10 8 s -1 .
The photodegradation rate k d is set to 170 s -1 . Using these values it is possible to calculate the parameters that enter into Eq. 31, b = σ P c /k r x = 3.117 ⋅ 10 -6 cm 2 s/J and η = k d b = 5.3 ⋅ 10 -4 cm 2 /J. The solution of Eq. 31 is written as where L = 0.02. Equations (21) and (24) can be used to calculate α = 6324 and α′ = 8462, respectively, using Eq. (32). Note that A is a constant which is set to 1. Figure 3 shows the resulting N(x,y) for the parameters given above. The solution given by Eq. (32) can be used to predict the time independent fluorescence signal by substituting it into Eq. (19). The result can be compared to fluorescence calculated from the steady state of the approximate kinetic Eq. (30), i.e., The two values of the fluorescent signal (FS) agree to within 1 % up to 0.1 W incident power. The approximate fluorescence signal given by Eq. (30) and Eq. (33) is almost indistinguishable from the result given by Eq. (32). It is expected that the time dependent solutions will also be described correctly by the approximate Eq. (30) which together with Eq. (19) will be taken as the fundamental description of the photodegradation process occurring in the focused laser beam. Finally the steady state solution was used to find the effective distance to use in the estimate of the average spatial derivative given by Eq. For the case of low power, the approximations bα P(t) « 1, bα′ P(t) « 1, and a laser beam power given by P(t) = P 0 + P 1 cos (ω t) yield an approximate solution given by Eq.    N(x,y). The derivative is finite along the profile of the laser beam.  Figure 5 shows the time dependent solution at three values of the laser power for a modulation frequency of 2 Hz. The other parameters were taken to be the same as in the case of the steady state solution. The salient features of the time dependent solution are a transient at t = 0 and a modulation of the intact fluorophore concentration at the driving frequency of the laser. As expected the average level of intact fluorophores decreases with increasing average power, and the amplitude of the modulation increases with increasing modulation amplitude. The approximate time dependent solution given by Eq. (36) was used to obtain the fluorescence signal given by Eq. (19). The fluorescence signal was detected by a computer model of a lock-in amplifier consisting of an in-phase and quadrature multipliers followed by low pass filters. The low pass filters were implemented by Fourier analyzing the signal, setting to zero all coefficients corresponding to non zero frequencies, and then performing an inverse Fourier transform. The solid circles in Fig. 6a show the ratio of the detected quadra-ture and in-phase components at different modulation frequencies. Most of the parameters used in Eq. (36) were the same as in the steady state calculation except for width = 10 µm, P 0 = 50 mW, and k D = 600 s -1 . The solid line in Fig. 6a shows a fit to the calculated points of a function given by R( f ) = (a* f )/(b + f 2 ), where f is the modulation frequency and a and b are fit parameters. Clearly the calculated response and R(f ) form an excellent match. The solid points in Fig. 6b show experimentally measured values of the ratio of quadrature and inphase response. The solid line in Fig 6b is a fit to the function R(f ). The match is reasonable. Thus there is a strong indication that the measured response can be represented by the function R(f ), and that the fit parameters provide information about the various physical constants of the experiment. To elucidate the dependence of the fit parameters on the physical constants, it is best to use analytical solutions to Eq. (30) and obtain the analytical form of the response function R(f ). This task is carried out in a forthcoming paper. We close this section with a brief description of the mathematical justification for the time dependent case. The reduction of the original system Eq. (4) and Eq. (5) to a single partial differential equation depended on the use of a multiple time scale perturbation expansion. The validity of the technique used for this problem (matched asymptotic expansions) is discussed in Appendix A. By spatially averaging the partial differential equation we derive an ordinary differential equation for the averaged reactant concentration as well as a simplified expression for the total fluorescence. These simplifications are possible because the varying functions (Eq. (21) and Eq. (24) respectively) are nearly constant. The proof of this fact can be found in Appendix B. Finally the equation for the spatial average, as in Eq. (30), involves these functions and it is shown that replacing them by constants creates a small error. This is shown in Appendix C.

Conclusion
The dynamics of the populations of excited and unexcited fluorophores are described by a pair of first order partial differential equations that incorporate the kinetics of the fluorophore reaction driven by a periodically varying laser beam and the convective fluid flow in a measurement apparatus. We obtain a dramatic reduction in the mathematical complexity of the model by taking advantage of the disparate time scales for excitation, relaxation and degradation. This is followed by averaging the reactant concentration over the laser beam, which is assumed to have a Gaussian profile. A single ordinary differential equation, Eq. (30), is obtained describing the evolution of the spatially averaged reactant concentration. The measured fluorescent signal, Eq. (19), is expressed in terms of the spatially averaged reactant concentrations. These equations provide a prediction of the response whenever the photochemical reaction degrades the fluorophores flowing past the laser beam.
Mathematical justification for the approximations used in the course of the derivation is given and some comparisons to the time independent case are described. The predicted frequency response (based on the time varying solution of Eq. (30)) is used to derive a functional form for the ratio of the in-phase and quadrature parts of the fluorescent signal. This same form provides a good fit for experimental data. The fit and the connection with physical parameters of the experiment will be discussed in a forthcoming paper.

APPENDIX A
This appendix addresses the issue of the validity of the perturbation expansion of N 0 and N 1 for large times as stated in Eq. (6) and Eq. (7). The discussion is based on the standard observation that the solutions of Eq. (6) and Eq. (7) can be expressed in terms of the solution of a parameterized system of ordinary differential equations. Here a new time variable s where t = s + θ is introduced. The relaxation rate is denoted by k r x .

(A1)
See for example the book Aris and Amundson [9], which contains many worked examples from the chemistry literature as well as a general discussion of the transformation of partial differential equations such as Eq. (6) . This is known as the "method of characteristics." The variable θ ≥ 0 parameterizes the family of solution curves that satisfy the differential equations in s and the initial conditions at s = 0 displayed on the last line. The original time variable is t = θ + s and this leads to the interpretation of θ as the time a particular family of particles entered the apparatus (at x = -L / 2 ) while s is the amount of time elapsed since entry. If we are following the set of particles that entered the apparatus at the beginning of the observation period, i.e., at θ = 0, then s = t. Here we will assume that θ ≥ -L / v . For each θ, a perturbation method will be used to construct The appropriate initial conditions can be obtained by examining the behavior of the inner corrections. Equations for them are obtained by returning to A1 and replacing the functions N 0 , N 1 , and x and by the expressions on the right hand side of A2. We seek representations of the corrections in terms of the following perturbation expansion: We will suppress the dependence on y in what follows. After performing the previously described replacement involving A1and A2, and then using the right hand sides of A4, the desired equations are obtained. In them, the fast time variable is introduced by the change of variable, s = ετ. We will write down the resulting equations just for η 1 and ξ, The equation for η 0 is quite similar to the one for η 1 .
The initial conditions of the original problem A1 can be used as the initial conditions for η 0 ,η 1 , and ξ since these functions approximate the initial part of the solution of A1. Thus The initial value X 0 (0) is a constant that will be determined shortly. Adding the first two equations in terms in the expansion of the inner correction of N 0 + N 1 is a constant. The next step in determining the initial conditions for the large time solution (which in general will not be the initial conditions for the original problem) is to match the large time or outer solution with the inner correction in an overlapping region of the s axis. This region is defined as the outer reaches of the initial interval where the inner correction is valid and the beginning of the region of validity of the large time solution. To do this we calculate the limit of the inner correction as τ → ∞ by fixing s and letting ε → 0. The limit must match the limit of the large time solution as s = ετ → 0 where τ is fixed and ε → 0. For our problem this is written as, In this way, we determine the initial condition for the leading order large time solution N(y,θ,s) and thereby determine the boundary condition. Indeed from the second equation in A10, the initial condition for the x position to leading order is X 0 (0) = -L / 2 . Given the large time solution and inner corrections formally constructed here, we can for each θ ≥ 0 construct a solution of A1 that is uniformly valid throughout the interval 0 ≤ s ≤ T /ε for each fixed ε > 0 Standard theorems in the theory of singular perturbations show that the operations used here to obtain it are rigorously justified. The key hypothesis that must be satisfied in our problem is that the inner solution must approach a finite limit exponentially. Recall that η 10 = N 0η 00 . Substituting this relation into the first equation in A9 results in the following equation for η 00 . (A11) Solving this constant coefficient equation shows that η 00 and thus η 10 has the required behavior as τ → ∞. Note that equation (A1) is not a standard singular perturbation problem because the Jacobian associated with the reduced problem, i.e., the problem obtained from setting ε = 0 in the first two equations of (A1), is singular. Validity of the singular perturbation method in this "critical" case was proved by Vasil'eva and Butuzov [10] and was treated analytically and numerically by O'Malley and Flaherty [11]. These results require that the Jacobian have a constant deficient rank k for all s and its non-zero eigenvalues have negative real parts. A discussion of these problems can be found in O'Malley [12] and Smith [13].
A complete solution of (A1) uniformly valid in an interval comes from combining the large time and inner correction and then subtracting the common value they share in the overlapping region. If N(y,θ ,s) is the total fluorophore population to leading order we have, where δ (s) and β (s) are exponentially small in s and O(ε) holds uniformly in the interval 0 ≤ s ≤ T /ε and where X 0 (s) = -L / 2 + νs. On adding Eqs. (A13) and (A14) we see that at large times, the total fluorophore population is dominated by N(y,θ,s). It is in this sense that we can justify the statements in Sec. 2.  We assume that all of the limiting averaging operations used to derive B1 are valid. This will be the case for example when the integrands are quasi-periodic  The conditions for Ψ at x = -L / 2 , are obtained from the conditions for Φ and N. Note the discontinuity in the value of Ψ at t = 0 as particles are pumped into the vessel. It propagates to larger values of x for t > 0 as a shock. However since here we will assume that will have passed and we have continuity in Ψ. The (B3) respect to x is written, where for easier reading we have suppressed some of the function arguments. The functions H t and g t are derivatives (B5) From Eq. (B2) it can be seen that g t = -H t | ξ = x′ Φ(r).
show that g, H, H τ are all small. Indeed from Eq. (B1) we have Φ(r) ≤ N 0 so | g t | will be small if | H τ | is. Using the fact that P(r,t) = (P 0 + ∆(t)) f (r) where ∆(t) = P 1 cos(ωt) we find that (B6) We can perform a similar analysis of the remaining terms in Eq. (B3) and Eq. (B4). Using the notation O(ρ) to denote a term whose absolute value remains bounded by a multiple of ρ as ρ → 0, the results are summarized as