Numerical study of a convective cooling strategy for increasing safe power levels in two-photon brain imaging

Two-photon excitation fluorescence microscopy has become an effective tool for tracking neural activity in the brain at high resolutions thanks to its intrinsic optical sectioning and deep penetration capabilities. However, advanced two-photon microscopy modalities enabling high-speed and/or deep-tissue imaging necessitate high average laser powers, thus increasing the susceptibility of tissue heating due to out-of-focus absorption. Despite cooling the cranial window by maintaining the objective at a fixed temperature, average laser powers exceeding 100–200 mW have been shown to exhibit the potential for altering physiological responses of the brain. This paper proposes an enhanced cooling technique for inducing a laminar flow to the objective immersion layer while implementing duty cycles. Through a numerical study, we analyze the efficacy of heat dissipation of the proposed method and compare it with that of the conventional, fixed-temperature objective cooling technique. The results show that improved cooling could be achieved by choosing appropriate flow rates and physiologically relevant immersion cooling temperatures, potentially increasing safe laser power levels by up to three times (3×). The proposed active cooling method can provide an opportunity for faster scan speeds and enhanced signals in nonlinear deep brain imaging.


Introduction
Recent advances in two-photon excitation fluorescence (TPEF) microscopy [1] have enabled deeper and faster neuronal imaging to better understand how the brain processes information and reacts to stimuli [2][3][4][5][6].The inherent optical sectioning capabilities [7][8][9][10] make TPEF especially suitable for in vivo brain imaging.Specifically, superficial structures that are accessible in vivo with minimal craniological invasion, such as the neocortex and cerebellum [11], are often of interest to TPEF microscopy.Since the probability of simultaneous absorption of two photons is low, high intensities become desirable.On the other hand, high-speed imaging deep in tissue requires high average powers, which make the tissue susceptible to thermal damage.While bulk heating depends on average power [12][13][14], high peak intensity pulses may elicit nonlinear photodamage [15][16][17].Beam focusing introduces the chances of thermal accumulation during scanning not only within the focal volume but also in the out-of-focus region.More importantly, the light intensity reaching the focus drops exponentially due to tissue scattering, requiring an increase in average powers with depth to maintain desired signal levels.Such an increment in power increases the chances of out-of-focus absorption and bulk heating.High-speed and high-resolution deep brain imaging is, thus, limited by the risk of tissue heating and thermal damage.
In modern neuroscientific applications, where two-photon imaging is used as a tool to track neuronal action potentials and calcium dynamics through an entire volume, the problem of tissue heating has yet another facet.Past studies have demonstrated changes in neural activity and spikes in action potential due to thermal gradients induced by photon absorption [12,[18][19][20][21].Such changes may cause artifacts in the brain activity to be recorded.As a result, the characterization of heating and thermal transport becomes important in the realm of two-photon imaging [12,13,[22][23][24][25].
The heat dissipation problem in brain imaging is non-trivial with a limited scope of existing cooling mechanisms.Till date, neuronal imaging experiments with a cranial window have only relied on passive cooling via a cooled objective in conjunction with or without active cooling via blood perfusion mechanisms for heat dissipation [13,26].While this conventional method of cooling is able to alleviate thermal accumulation in the tissue to some degree, laser powers have still been limited to 100-200 mW.This is due to the passive nature of maintaining the objective at a fixed cold temperature, where heat dissipation saturates over time.External cooling also makes the superficial layers of the tissue susceptible to a substantial decrease in temperatures.Although, it is known that the brain is neuroprotective in nature to hypothermia, recovery from changes in blood flow and oxygenation due to temperature drop takes place over a span of days [27].
Other cooling approaches have relied on reducing the laser repetition rate.While brain imaging depths of up to 1 mm were successfully reached [28], there are consequential limitations, in general, with reducing the repetition rate [29][30][31].Among others, repetition rate reduction imposes a limitation on imaging speed, which is especially a cause of concern in fast volumetric imaging required to capture rapid brain activity [32].Thus, there is a clear need to devise enhanced cooling strategies supported by a comprehensive evaluation of thermal transport.
Here, we study an alternative cooling method that relies on convective heat transfer by inducing fluid flow, which fundamentally poses the potential of further enhancing thermal diffusion within the tissue.The objective of this study is to theoretically propose and investigate this active cooling technique for deep brain imaging in a craniotomy setting to enable higher speed two-photon microscopy.A comprehensive parametric study compares the proposed convective cooling method with induced laminar flow to the conventional method with immersion water remaining stagnant.A parameter space consisting of laser excitation powers, imaging depth, numerical aperture (NA), and field of view (FOV) has been swept for presenting an in-depth analysis of the effects of the said parameters on thermal accumulation for an excitation wavelength of 1035 nm in a point scan modality.Furthermore, we are presenting a detailed account of the impact of applying duty cycles from a thermal dissipation point of view.Numerical simulations show that the proposed convective cooling leads to an increase of average excitation laser powers by 3× over the currently established safe power limits with conventional cooling methods.Comparisons with studies available in the literature indicate that this cooling method has the potential to improve the brain imaging depth by 50% at 920 nm wavelength without compromising imaging speed while maintaining physiological temperatures.The simulation software was written in MATLAB (Mathworks Inc.).

Methods
The physics of light transport within a biological medium can be examined with the help of the Radiative Transport Equation (RTE).The RTE is an integrodifferential equation [33], which is complex to solve for the case of a participating medium.Analytical solutions such as the diffusion approximation can provide computationally efficient solutions [34].However, the diffusion approximation is not applicable for moderate to highly scattering media.Instead, the present study uses the Monte Carlo (MC) technique that provides a gold standard pathway for numerically investigating optical transport through tissue.
We adapted the rules for performing the MC simulations from the MCML package [35] with an improved approach for Gaussian beam propagation, inspired by a spot focusing technique [26,36].We sampled photons on the tissue surface in a Gaussian distribution based on the 1/e 2 beam diameter, objective NA, and imaging depth.Trigonometric formulations determine each of the photons' initial coordinates and propagation directions right before their launch.Launch coordinates (subscript C) and directionality (subscript D) of each of the photons are given by Eq. ( 1) and Eq.(2). (1) (2) Here, r i is the initial photon radial position at the surface, f is the objective focal length, n is the tissue refractive index and ξ is a pseudorandom number.θ and ∅ are phase and azimuthal angles, respectively.The photons take pseudorandom finite step sizes that are derived from an inverse transform applied to the Beer's law.An event of interaction takes place only when the photon has traveled by a distance, s ξ , tantamount to the predetermined step size based on the wavelength-dependent total attenuation length (Eq.( 4)).
Here, µ t is the wavelength-dependent total attenuation coefficient.At the site of interaction, the photon weight is updated by subtracting ∆W from its current weight, owing to absorption in the tissue medium (Eq.( 5)).This subtracted weight is added to the weight of the corresponding pixel.
Here, W is the current weight of the photon and µ a is the wavelength-dependent tissue absorption coefficient.After the weight update step, the photon changes its direction of trajectory based on the Henyey-Greenstein scattering phase function (Eq.( 6)).
Here, g is the tissue anisotropy coefficient.A game of roulette checks if the photon survives at the interaction site.Figure 1(a) presents a flowchart for explaining the general workflow of the MC simulations.Whenever a photon reaches a location that is at or beyond the tissue domain boundaries, we relocate the photon at the surface without any weight.Further mathematical details of the rules followed for carrying out the MC simulations and the spot focusing technique may be found in the works of Wang, Jacques, and Zheng [35] and Blanca and Saloma [36].Implementing a time-dependent MC model for simulating FOV scan is computationally expensive.Past experimental findings [13] demonstrate that the brain tissue temperature varies only by 0.1-0.15°C/sunder 20 s of excitation for a FOV scan of 1 mm, post which the temperature reaches steady state.For high frame rates, therefore, such a temperature variation permits modeling of the FOV scan as a steady-state phenomenon simplifying the MC simulations greatly.In our simulations, the beam profile is Gaussian at the time of launch and the Gaussian beam profile is randomly distributed across the FOV scan area on the surface.The amount of heat accumulation depends on various laser dosimetry parameters, FOV size, and the NA.Blood flow within the arteries and veins acts as a convective cooling source within the tissue medium besides the presence of the metabolic heat generation source.Essentially, we obtain the fluence rate distribution using the MC technique and use it as a source term in the Pennes' bioheat transfer equation [37], given by: where, ρ, c p , k, and T are the density, specific heat, thermal conductivity, and temperature field of the tissue, respectively.t and (x, z) represent time and spatial coordinates, respectively.Subscript bl denotes the thermophysical properties of blood and T A is the arterial blood temperature.ω bl is the blood perfusion rate.q m and ϕµ a are the metabolic heat generation and laser fluence rate intensity source terms, respectively.We discretize the bioheat transfer equation in space and march it in time using the explicit finite difference method (FDM), specifically by using central differences for the first-and second-order spatial derivatives and forward differences for the time derivative.We performed our heat transfer simulations on a computational domain consisting of the brain tissue, skull bone, cover glass, and water immersion layer, as shown in Fig. 1(b), which also indicates the relevant boundary and initial conditions.The objective front aperture diameter is greater than 8 mm.Thus, we applied a cold thermal boundary condition on the entirety of the water medium.In a physics-driven numerical simulation such as the one currently under consideration, the accuracy of the predicted results largely depends on the grid discretization.As a result, it is imperative to perform a reliable grid dependence analysis, before selecting the spatial and temporal mesh sizes.To determine the appropriate grid size, we continuously reduced spatial and temporal step sizes until changes in the predicted results ceased to exist.We found the grid independent mesh size and temporal step size to be 0.0285 mm and 0.15 ms, respectively.Relevant thermophysical properties were taken from [38,39].
We considered two different cooling strategies.First, for the simulation of conventional objective cooling, where the objective is held at fixed temperature, we maintained a constant boundary temperature (Dirichlet boundary condition) at the top surface of the water medium.For the glass surface at the interface with water, we applied a linear distribution of temperatures ranging between the water and body core temperatures.The glass surface at the interface with tissue has a constant initial temperature of body core and we assumed a linear distribution between the surfaces.Second, we investigated the proposed convective cooling method by inducing a laminar flow to the immersion water layer between the objective and the cranial window.Such a scheme requires the application of the Robin boundary condition (Eq.( 8)), which relates convective and diffusive transport across an interface.For implementing this cooling scheme, we maintained the entire water medium at a fixed temperature and applied the following boundary condition at the water-glass and water-bone interfaces: Here, h and k are the convective and conductive heat transfer coefficients, respectively, and T e is the flowing water temperature.However, one caveat to note with regards to the flowing water-cooling scheme is that in addition to temperature being the scalar to solve for, velocity and pressure terms also get introduced into the problem.This may be explained as follows.The physics of flowing water is governed by the continuity and momentum equations with pressure and velocity being the unknown scalars.Solution of the pressure and velocity terms in the continuity and Navier-Stokes momentum equations requires a suitable technique to decouple the pressure and velocity terms.To determine the temperature field at the interface, the velocity field obtained from the said numerical procedure would then need to be applied to the general transport equation for temperature using: Here, v x and v z represent the velocity field and S g is the source term.While the mentioned solution procedure is an ideal candidate from a theoretical perspective, the computational power and time required to perform the same would be significant due to the additional equations that need to be solved.As a result, we made an assumption from analytical relationships that simplified the problem.We determined the permissible flow velocities from past studies involving water irrigation (achievable flow rate of Q = 10-15 mL/min) in laser ablation applications [40] and constrained the Reynolds number (Re) within the laminar regime.Based on the Nusselt number (Nu = 7.54) for internal flows with high aspect ratio, we obtained a range of suitable convective heat transfer coefficients (h = 1-5 mW/mm 2 -K) [41].These values of convective heat transfer coefficients were applied at the Robin interface to enhance heat diffusion within the tissue, as would be seen in the upcoming section.
Here, ρ and µ are density and dynamic viscosity of water, respectively, and D and A are the characteristic length and area of the water contact region.The rationale behind applying such an assumption is that the time scale required for heat to penetrate into the immersion water layer is much larger than the time it takes for water to flow away.As a result, the water layer may always be statically thought of as being a cold medium with convective thermal dissipation occurring only at the interface.
Another simplification that can be applied to heat transfer simulations with high repetition rate pulsed sources is that the laser excitation can be assumed to be of continuous wave type if the thermal relaxation time (τ R ) of the tissue, given by Eq. ( 12) [42], is much larger than the time between two pulses.
Here, d is the beam spot diameter and α is the tissue thermal diffusivity.In the case of brain tissue, the heat diffusion time scale (µs) is three orders of magnitude larger than the time between pulses (ns), which renders the problem as being that of continuous heat buildup.The stated simplification thus holds valid in our case setup.

Results
This section presents the findings of our parametric study to compare the efficacy of the proposed convective cooling and the conventional objective cooling techniques over a parameter space consisting of imaging depth, average surface power, water temperature, FOV, and NA.We also evaluate the effect of implementing duty cycles and present an in-depth analysis of the proposed cooling technique combined with the repetition rate reduction strategy for facilitating brain imaging at enhanced speeds.Our simulations consider an imaging setup with an Ytterbium fiber laser (Monaco, Coherent Inc.) as the excitation source at 1035 nm.Such laser systems are particularly suitable for high-speed imaging [4,43,44] as they can provide tunable repetition rates ranging between 1 Hz and 50 MHz at the desired high average powers.From here on, we shall refer to the conventional cooling technique as Dirichlet cooling and the proposed convective cooling as Robin cooling.

Model validation studies
For beam focusing onto the sample, we considered a water dipping, high NA objective (Olympus XLUMPFL 0.95 NA, 20× magnification) having a working distance of 2 mm.The 1/e 2 beam diameter at the objective back aperture was 12 mm, giving a fill factor of 70% [45].In order to verify our mathematical model, we compared our numerical results with thermocouple contact thermometry-based experimental brain heating results from the work of Podgorski and Ranganathan [13].The excitation wavelength and cooling technique along with all other aforementioned optical parameters were taken from their work.We found an error of less than 10% between our simulations and the experimental data in terms of maximum temperature change per 100 mW at both 150 µm and 500 µm imaging depths for all of the three wavelengths studied, viz., 800 nm, 920 nm, and 1040 nm, which gave our simulations a high confidence value.Figure 2 summarizes the results of our model validation studies.Fig. 2. Numerical verification of our model in terms of peak temperature rise per 100 mW of average laser power at 150 µm and 500 µm imaging depths for 1040 nm, 920 nm and 840 nm wavelengths.Experimental data was obtained from the work of Podgorski and Ranganathan [13].E: Experimental data from [13], S: Simulation results.

Effect of convective cooling
We compared the two cooling techniques for six values of cold boundary temperatures (immersion water temperature in Dirichlet case and flowing water temperature in Robin case) between 23°C and 34°C.Robin cooling with the help of water recirculation is expected to provide enhanced thermal takeoff near the surface when compared to Dirichlet cooling.Figure 3 shows this degree of enhancement for a FOV of 500 µm and NA of 0.95 with a representative boundary temperature of 25°C.In these plots, we represent imaging depth in terms of wavelength-dependent total attenuation length, l t (0.168 mm for 1035 nm in gray brain matter [46]).The variation of temperature (after the objective) with excitation power exhibits a linear dependence of maximum temperature on laser excitation power (Fig. 3(a)) in agreement with previous studies [13,26].Although high average surface power values are not required for imaging at the surface and superficial regions, we plot the entire parameter space of imaging depth and powers in Fig. 3(a) to reflect laser powers that might be needed for ultra-high frame rates (> 10s of kHz).
The peak temperatures seem to be generally higher when focusing near the surface than at deeper extents for the same power in case of Dirichlet cooling.However, some exceptions evolve in case of Robin cooling, as may be noted from the temperature curve corresponding to 0l t .The reason may primarily be attributed to two factors.The first is the proximity of imaging depth to the cooling front, which is favorable for better heat diffusion.In Dirichlet cooling, the distance between the cooling front (objective front aperture) and the focal plane does not change with imaging depth.Consequently, the proximity factor remains constant for all the imaging depths.In Robin cooling, on the contrary, the interface between flowing water and tissue surface acts as the cooling front.The distance between this interface and the focal plane keeps increasing with imaging depth.The second factor concerns the deposited fluence.With increasing imaging depth, there is better delocalization of fluence.Combined with superior thermal take off at the convective interface with Robin cooling, this delocalization factor leads to enhanced overall heat dissipation from in and around the focal volume.Thus, Robin cooling facilitates lower peak temperatures for both small and large imaging depths in comparison to moderate imaging depths.
To quantify the thermal effects induced on the tissue during imaging, we defined a heat affected zone (HAZ) parameter.The physiological temperature range of the brain has been defined by the threshold temperature values of 34°C and 39°C inflicting hypo-and hyperthermia, respectively [18,47].The HAZ, therefore, represents the total summation of annulus cylindrical pixel volumes where temperatures exceed 39°C.The variation of HAZ with power shows similar trends to those of peak temperatures.The proposed Robin cooling technique indeed provides superior heat dissipation than Dirichlet cooling, as can be seen by the lower peak temperatures in Fig. 3(a) and smaller HAZs in Fig. 3(b).It would be most ideal to be able to perform imaging without the occurrence of any HAZ.The HAZ plots show that convective cooling with h = 5 mW/mm 2 -K doubles the power allowance (for zero HAZ) for a given imaging depth.On the other hand, Robin cooling facilitates up to 50% reduction in HAZ for a given power value.The trends remain similar irrespective of the cold boundary temperature (data not shown here).Thus, Robin cooling is seen to greatly enhance heat dissipation in and around the focal volume.
Interestingly, Robin cooling subjects the sub-surface regions to significant cooling due to better proximity to the cold interface such that temperatures drop below the previously mentioned physiological range.We further characterize cooling schemes in terms of cooling affected zone (CAZ) by calculating the total volume where temperature drops below 34°C at the end of 60 s of imaging.Variation of CAZ with power shows that there is a larger CAZ associated with convective cooling (Fig. 3(c)).Due to increased thermal accumulation at higher powers, the CAZ is seen to decrease with excitation power.Therefore, while Robin cooling poses the potential of using enhanced imaging powers from the perspective of reduced HAZ, the risk of a larger associated CAZ calls for more consideration in averting excessive sub-surface cooling.
Table 1 presents the maximum average power values that ensure zero HAZ (P max, HAZ=0 ) for various water flow temperatures, T e , at 6l t imaging depth.We will see in section 3.4 that these power values could be increased with implementation of duty cycles.We also show the CAZ values so that scientists can make an educated choice on flowing water temperature depending on their specific experimental conditions.Figure 3(d-g) show representative temperature contours for the two cooling schemes at the end of 60 s of laser excitation for two imaging depths viz., 0l t (surface) and 6l t at 50 mW and 500 mW power, respectively.It is immediately evident from the temperature distributions that the efficacy of thermal dissipation is much better in the case of Robin cooling than Dirichlet cooling.The lateral extent of the region having potential for undesirable thermal effects is greater in Dirichlet cooling due to heating of the stagnant immersion water layer.Previous numerical studies incorporating cooling with a water immersion objective held at a fixed temperature have reported a shift of thermal accumulation region from the focal volume [13].We intended to compare such a spatial response for the Robin cooling scheme against the Dirichlet cooling technique.Figure 3(h) shows variations of temperature with depth along the central axis of the tissue domain at 6l t with a cold boundary temperature of 25°C for both the cooling schemes.For reference, the location that corresponds to the focal depth has also been marked on the plots.It may be observed that the heat accumulation region is above the focal plane for large imaging depths.This result is primarily attributed to the factor of increased absorption above the focal plane at high imaging depths.On the other hand, for shallow imaging depths, we observed a tendency of the heat front to diffuse downwards with respect to the focal plane (data not shown here).
Figure 3(i) shows the temporal variation of temperature at the center of the FOV for 6l t .Robin cooling with high convective heat transfer coefficient (h = 5 mW/mm 2 -K) takes the least time to stabilize the temperature at the focus.This effect is especially pronounced at shallow imaging depths.An important point to note is that the disparity between Dirichlet cooling and Robin scheme with h = 1 mW/mm 2 -K in terms of spatial and transient thermal response almost closes out with increasing imaging depth.Thus, there is little advantage in using h = 1 mW/mm 2 -K, especially at large imaging depths.
It is important to note that the cold boundary temperature plays a major role in CAZ formation [27].The intermediary period between finalization of the experimental setup for imaging and the start of imaging could play a detrimental role by causing the sub-surface tissue region to approach the cold boundary temperature.This delay, in fact, plays a major role in creating the large CAZs that were seen in Fig. 3(c).An approach to tackle the said problem could be to use a cold boundary temperature that is higher than the room temperature and corresponding to the lower bound of the physiological range, which would lead to substantial reduction or even elimination of the CAZ.Such a discussion will follow in the upcoming sections, where we consider cold boundary conditions ≥ 32°C.

Effect of FOV and NA
The objective of this section is to compare responses of both cooling schemes for various FOVs and NAs.As it would intuitively be expected, depositing the same amount of power over a large FOV would lead to less heat accumulation per unit area.Previous reports have shown that it is more convenient to image a large FOV from the heating perspective [26] due to enhanced thermal transport by diffusion.Our simulations with the Dirichlet cooling scheme at all of the cold boundary temperatures and imaging depths, verify this finding.Figure 4(a) and 4(b) compare the temperature distributions of the Robin scheme with those of the Dirichlet scheme for a cold boundary temperature of 32°C at two imaging depths: at the surface (0l t ) and at 6l t using 50 mW and 350 mW average powers, respectively.The Robin scheme demonstrates a similar spatial response with FOV variation when compared to the Dirichlet scheme in that, peak temperatures near the focal volume decrease with an increase in FOV.We found a scaling relationship (Eq.( 13)) between peak temperatures (T peak ) for a given FOV and average power.The entire parametric dataset consisting of imaging depths (z), average powers (P s ), flowing water temperatures, FOV, and NA was studied for the Robin cooling scheme to determine this empirical relationship.
Figure 4(c) presents the effect of varying the focusing NA on tissue heating at an FOV of 1 mm for the Robin scheme at a flowing water temperature of 32°C.It is important to note here that the NA of the system was reduced from 0.95 to 0.6 in our simulations by underfilling the objective while keeping the working distance unchanged.Results show that reducing the NA brings down the average power ceiling for ensuring zero HAZ by up to 25%.This is intuitive because reducing the NA at the same imaging depth decreases the focusing cone angle within the tissue.The increment in photon density within the focusing cone creates more chances for absorption and a stronger fluence rate source term as a consequence.These results indicate that there is less susceptibility to heating with imaging at high NA.

Effect of duty cycles in Robin cooling scheme
Implementing laser ON-OFF duty cycles during imaging is another potential pathway towards alleviating thermal accumulation within the tissue [13].The rationale behind such an approach is to allow some time for heat to diffuse away between each laser excitation cycle and prevent accumulation.The characteristic thermal diffusion time for the brain tissue is on the order of 1 µs.A laser OFF time greater than the mentioned characteristic time scale would allow for thermal transport through diffusion.Duty cycles can be implemented in an imaging setup with the use of programmable mechanical beam shutters.The cycle timing is therefore limited, for example, by the mechanical lag time of the shutter.The fastest mechanical shutters have a response time of 5-10 ms, which is another factor to consider when designing the duty cycle.To demonstrate the effect of duty cycle on heat accumulation for the Robin cooling scheme, we performed our simulations for a case of 0.5 mm FOV and 0.95 NA at 50% duty cycle.Here, 50% pertains to the laser ON time.We specifically simulated two cycle times of 0.1 s and 1 s, where 50% duty cycle would refer to a laser ON time of 0.05 s and 0.5 s and a laser OFF time of 0.05 s and 0.5 s, respectively.
Figure 5(a) depicts the effect of duty cycle on heating at a cooling boundary temperature of 32°C.The advantage of implementing duty cycles can be seen from the variation of peak temperature within the tissue domain.Changes in peak temperature with respect to body core temperature at large imaging depths drop by up to 70% with duty cycles.These temperatures rise above 39°C at 500 mW when a cycle time of 0.1 s is used.On the contrary, this power threshold drops to 400 mW with a cycle time of 1 s.As an additional observation, the disparity between the peak temperature curves corresponding to various imaging depths reduces significantly when duty cycles are implemented.
The transient response of temperature as a result of duty cycle shows an intermediary diffusion phenomenon occurring between laser excitation cycles and an overall reduction in temperature at the focal plane (Fig. 5(b) and 5(c)).The same trend is true for any location within the tissue domain.However, the temperature fluctuates within ±1.5°C at 50% duty cycle with 1 s cycle time.Reducing the cycle time can effectively minimize or prevent such temperature rise and fall patterns, as we may observe from the transient response curves with fluctuations minimized to ±0.5°C for 0.1 s cycle time (Fig. 5(c)).
The HAZ and CAZ plots provide a more holistic picture of the overall effects of implementing duty cycles (Fig. 5(d)).The power up to which HAZ is seen to remain non-existent becomes more than two times (2×) with 50% duty cycle at a 0.1 s cycle time for all imaging depths.Furthermore, using a cold boundary temperature of 34°C (lower limit of the physiological range) enables the elimination of CAZ all together.The HAZ plot shown in Fig. 5(d) indicates that a maximum of 350 mW average power can be applied at the surface for imaging 1 mm deep (6l t ) at 1035 nm with the optical parameters under consideration.

Outcomes of enhanced cooling and discussion
To quantify the advantages of implementing Robin cooling alongside duty cycles for two-photon deep brain imaging, we present a comparison between thermal aspects resulting from conventional Dirichlet cooling and those from Robin cooling (from hereon, referred to as the proposed cooling technique) with the application of 50% duty cycle of 0.1 s cycle time to both the schemes.This study considers optical parameters taken from the work of Wang et al. [26] (1.05 NA with ∼75% back-aperture filling, effective NA of 0.75, 7.2 mm focal length, 2 mm working distance, 0.3 mm FOV, 920 nm wavelength, l t = 148 µm, 60 fs pulse width, 80 MHz nominal repetition rate).For these optical parameters, a pulse energy of 0.2 nJ is required at the focus to obtain two-photon fluorescence signal strength of 0.1 photon per pulse for GCaMP6s as the indicator [26].Due to scattering, the surface power needs to increase with imaging depth, potentially inflicting bulk heating effects in the sub-surface and focal region.An alternative strategy of varying the laser repetition rate can be used to maintain low average powers while increasing the pulse energy at the surface.In order to image ∼1 mm deep, the laser repetition rate is typically dropped to ∼1 MHz to avoid tissue thermal damage [26,28,48].The optimum laser repetition rate (f opt ) can then be estimated for safe levels of average power (P safe ) that limits thermal effects, and desired pulse energy at the focus (E f ), as a function of imaging depth (z) [49]: Our simulations and Wang et al.'s [26] results suggest a P safe value of 200 mW for the conventional cooling scheme beyond which, peak temperatures are seen to rise above 39°C.In contrast, our simulations suggest that we can achieve a 3× enhancement in average surface powers with the proposed cooling scheme, potentially allowing the safe delivery of 600 mW for deep brain imaging without the occurrence of HAZ and CAZ.Considering P safe = 200 mW as the safe average power limit, we plot the optimal repetition rate variation with depth in Fig. 6(a), while maintaining a constant pulse energy of E f = 0.2 nJ at the focus.Capped by a maximum laser repetition rate of 80 MHz (Tsunami, Spectra-Physics [26]) at the lower end of the imaging depth range, there is an exponential decrease of optimum repetition rate with depth.However, such reduction in repetition rate with depth may cause a reduction in imaging speed.Our analysis indicates that the proposed cooling technique with increased safe power levels enables imaging up to 550 µm depth, which is a 50% improvement, without any compromise in repetition rate and, thus, imaging speed.When compared to high-speed imaging without the application of any cooling mechanism (P safe = 100 mW [49]), we observe a 2× enhancement in imaging depth at uncompromised imaging speeds.
Without considering limitations in imaging speeds, the proposed cooling strategy with 3× enhancement in P safe enables an equivalent increment in E f for the same f opt variation as in conventional cooling.We evaluate the variation of average power after the objective, P s , based on imaging depth as per Beer's law and plot the same in Fig. 6(b).Scatter points indicate permissible power values at each imaging depth based on our simulations for both the cooling strategies for reference.The ability to deliver 3× the pulse energies could potentially enable a 9× enhancement in two-photon signal.Thermal contours presented in Fig. 6(c) for the case of our proposed cooling scheme with a flowing water temperature of 34°C demonstrate a visual proof of non-existent HAZ and CAZ.
Depending on laser dosimetry and beam focusing parameters, pulse energies might be potentially limited by photochemical damage thresholds.Comprehensive studies performed by Vogel et al. [42] indicate an irradiance value of 0.44 W/cm 2 for initiating photochemical damage.Based on the diffraction limited spot size corresponding to an effective NA of 0.75, we calculated the threshold pulse energy for initiating photochemical damage as ∼0.8 nJ.The correlation between initiation of photochemical damage and neuronal response remains to be understood since pulse energies of up to 1-2 nJ at the focus have been demonstrated to preserve neuronal physiological response under continuous stimuli [26].Pulse energies beyond this range may inflict ablation [50,51].
We note here that due to a smaller absorption coefficient at 920 nm, we were able to observe a 3× enhancement in average surface powers compared to 2× for 1035 nm wavelength.Notably, the theoretical power requirement is heavily dependent on laser pulse width besides other parameters highlighted above.As an example, if imaging were to be performed with 6 fs pulses [52] instead of 60 fs pulses, the power requirement would drop down by almost 70%.This may potentially provide further room for pulse energy enhancement limited by nonlinear effects at the focus and fluorescence saturation.

Discussion on implementation of enhanced cooling
In this section, we discuss some practical considerations on implementing the enhanced cooling strategy in craniotomy brain imaging experiments.Figure 7 presents a schematic of one possible example of incorporating a flow chamber above the cranial window.We envision a flow chamber that may be attached to the objective by a ball bearing mechanism, which allows the objective to be axially translated for imaging at different focal depths.Gaskets can be used for sealing of the objective.The other side of the flow chamber would have a window with a cover slip for the passage of laser excitation and emission.After creating a cranial window on the mouse skull, the animal head could be aligned with the window and fixed beneath the flow chamber, following common protocols in mouse brain imaging.The bottom surface of the flow chamber could also have a machined curvature based on standard radii of curvature of the mouse skull to minimize the distance between the coverslip and brain tissue.
While fluid structure interaction and associated thermal characteristics of the flow would be three-dimensional, and computational fluid dynamics simulations are recommended to predict them, one can make some two-dimensional analytical approximations considering the crosssection presented in Fig. 7.As an essential requirement to maintain steady flow characteristics in the region between the objective front aperture and coverslip, a fully developed flow over the window needs to be attained.Choosing a small value (say, 0.5-1 mm) of the hydrodynamic and thermal entrance lengths (L), we can calculate the Re value for flow over the coverslip, using Here, H img is the distance between the coverslip and objective front aperture.Pr is the Prandtl number and it may be taken as 6.9 for water within our temperature ranges.After determining the flow velocity, v img , from the calculated Re, we can determine the flow chamber height, H ch , by applying the Bernoulli's equation: Here, D obj is the physical outer diameter of the objective, g is the acceleration due to gravity (g = 9.81 m/s 2 ), and v ch is the flow velocity through the broader section of the chamber.In deriving this relationship, we considered the flow chamber width to be equal to the objective diameter, which may not be a necessary design condition.

Conclusion
The main goal of our study was to propose an active cooling technique of imparting laminar flow to the immersion water layer, facilitating a convective cooling pathway, for overcoming imaging speed limitations by increasing the safe power limits.We performed a rigorous parametric sweep at 1035 nm wavelength, where high average power ultrafast fiber lasers at high repetition rates are readily available to provide the necessary excitation powers for high-speed imaging deep in tissue.Through numerical simulations, we draw comparisons between the conventional cooling scheme of using a fixed temperature water immersion objective and our proposed technique of flowing the immersion water layer alongside implementing duty cycles in both the schemes.Our results indicate that convective takeoff at the interface of cooling water and glass coverslip facilitates enhanced heat diffusion within the tissue domain when compared to the conventional cooling technique.With duty cycle parametric studies, we show that the proposed cooling technique is especially effective for deep-brain imaging where there are high chances of out-of-focus absorption in the sub-surface region.We also discuss a possible implementation of the cooling strategy and present some relationships for determination of key design parameters of the flow chamber.
We quantified the possibility of causing an alteration in brain physiological response due to bulk heating and cooling by introducing the HAZ and CAZ metrics.To achieve minimal HAZ while avoiding undesirable consequences due to cooling, our studies show that the temperature of the circulation water should be kept close to the lower limit of the brain physiological range, which is, 34°C.FOV parametric studies indicate that imaging a large FOV allows for better heat diffusion.On the other hand, using a small focusing NA increases susceptibility to tissue heating.Lastly, we quantify the advantages of using convective cooling alongside duty cycles by performing a comparative study at 920 nm.The proposed cooling scheme poses the potential of enabling high-speed imaging with constant signal strength up to 550 µm imaging depth.Uncompromised imaging speed while maintaining good signal strength in brain imaging would especially be of relevance to ultra-fast imaging systems.On the other hand, our proposed cooling scheme enables a 3× enhancement in pulse energies at the focus and therefore, a 9× enhancement in achievable two-photon signal, when compared to the conventional cooling technique.We note here that certain dyes/probes used in functional brain imaging might have a narrower range of temperature sensitivity than the physiological temperature range defined in this work.Even in such cases, the key advantages of our proposed cooling technique would remain unaltered with the only condition that the flowing water temperature would be governed by the functional temperature range of the dyes/probes.
Besides two-photon imaging, the proposed scheme presents an opportunity of reducing the possibility of tissue thermal effects in three-photon imaging as well.Future studies involving both numerical and experimental characterization of thermal effects using the proposed cooling technique would help in better assessing its effectiveness on various other scanning and imaging modalities namely, line excitation and three-photon microscopy.

Fig. 1 .
Fig. 1.Simulation workflow and Monte Carlo results.(a) Flowchart demonstrating the workflow of the optothermal numerical solver.(b) Simulation domain with appropriate initial (subscript i) and boundary (subscript b) conditions.The domain size was 4 mm in each of the x and z directions.The glass slide radius and thickness were 2 mm and 140 µm, respectively.The bone thickness was 140 µm.The thickness of the immersion water layer changed with focal depth.For imaging at the surface, the thickness of the immersion was equal to the objective working distance of 2 mm.The idea of convective cooling was to impart a lateral flow to the water layer.(c)-(f) Normalized absorption profiles after a 0.2 mm FOV scan with 1035 nm at focusing depths of (c) 0 µm, (d) 336 µm, (e) 672 µm, and (f) 1,008 µm.

Fig. 3 .
Fig. 3. Comparison of Dirichlet and Robin cooling at 1035 nm using 0.95 NA and a scan FOV of 0.5 mm.Cold boundary temperature for both the cooling schemes was 25°C.(a) Variation of peak temperature with average surface laser power for Dirichlet scheme and Robin cooling at four different imaging depths.(b), (c) Variation of HAZ and CAZ with surface power for Dirichlet and Robin schemes with both h = 1 and 5 mW/mm 2 -K.(d), (e) Temperature contours after 60 s of laser excitation at the surface (0l t imaging depth) using 50 mW power for Dirichlet and Robin cooling schemes, respectively.(f), (g) Similar comparison of temperature contours for 6l t imaging depth using 500 mW power.(h) Temperature versus depth variation for 6l t imaging depth at 500 mW power.(i) Transient thermal response at the focal point for 6l t imaging depth at 500 mW power.All plots and contours presented in this figure show the thermal response at the end of 60 s of laser excitation.

Fig. 4 .
Fig. 4. Effect of varying the scan FOV and NA on thermal response with Dirichlet and Robin cooling.(a), (b) Variation of temperature with depth for the investigated FOVs at 0l t and 6l t using 50 mW and 350 mW of power, respectively.(c) HAZ formation for various NAs using Robin cooling with h = 5 mW/mm 2 -K and FOV = 1 mm.Cold boundary temperature of 32°C has been considered in these plots.

Fig. 5 .
Fig. 5. Effect of implementing duty cycles on thermal response with Robin cooling (h = 5 mW/mm 2 -K).The scan FOV was 0.5 mm with an objective NA of 0.95.(a) Comparison of no duty cycle and 50% duty cycle with 1 s and 0.1 s cycle times, respectively, in terms of peak temperature using a flowing water temperature of 32°C.(b), (c) Comparison of transient thermal responses at the center of the focal plane between cases of no duty cycle and 50% duty cycle with 1 s and 0.1 s cycle times, respectively.The flowing water temperature was 34°C and average powers of 100 mW and 350 mW were used for imaging depths of 1l t and 6l t , respectively.(d) Development of HAZ and CAZ with laser surface power.Arrows adjacent to the curves indicate relevant axes.CAZ values are zero throughout all the power values studied at each of the imaging depths.

Fig. 6 .
Fig. 6.Quantification of the benefits of using Robin cooling with duty cycles over the Dirichlet cooling scheme with parametric comparisons at 920 nm wavelength.The cold boundary temperature used was 34°C to ensure zero CAZ.Imaging parameters were taken from Wang et al. [26] for drawing comparisons.(a) Variation of optimum repetition rate with imaging depth for E f = 0.2 nJ at three P safe values of 100 mW, 200 mW and 600 mW.(b) Surface power variation with depth at the derived optimum repetition rate values (solid line plot in (a)) at P safe = 200 mW (E f = 0.2 nJ) and P safe = 600 mW (E f = 0.6 nJ).Permissible power values found from our numerical heat transfer simulations have also been overlaid for both conventional and proposed cooling schemes.(c) Temperature contours for the studied imaging depths are shown with no evidence of HAZ and CAZ.

Fig. 7 .
Fig. 7. Possible ideation of a flow chamber for implementing convective cooling.