Resolving the discrepancy between theoretical and experimental polarization of hyperpolarized 129 Xe using numerical simulations and in situ optical spectroscopy

For emerging biomedical applications of hyperpolarized xenon, the ability to obtain reliably high nuclear spin polarization levels is paramount. Yet, experimental nuclear spin polarization levels of xenon are highly variable and, more than often than not, well below what theory predicts. Despite rigorous and well-studied theoretical models for hyperpolarization and continuous-flow spin-exchange optical pumping (SEOP), there remains a substantial discrepancy between the theoretical and experimental polarization of 129 Xe; inexplicably, seemingly similar experimental parameters can yield very different polarization values. In this paper, the validity of the assumptions typically made about the thermodynamic state of the Rb vapor inside the optical pumping cell and the gas dynamics are investigated through finite element analysis simulations of realistic optical pumping cell models, while in situ optical and nuclear magnetic resonance spectroscopy measurements are used to validate the results of the simulations. Our results show that shorter xenon gas residence times and lower Rb vapor densities than those predicted by empirical saturated vapor pressure curves, along with incorrect SEOP parameters, are the primary cause of the discrepancy between theoretical and experimental polarization values reported in the literature


I. INTRODUCTION
Hyperpolarization of nuclear spins is used to increase the sensitivity of nuclear magnetic resonance (NMR) measurements by four to five orders of magnitude, thus enabling the detection of spin systems that otherwise would go undetected. 1 Among all the nuclear spins that can be polarized, 129 Xe presents unique features that make it an ideal probe for NMR investigations of materials and biological tissues, such as its inert nature and its high sensitivity to the chemical environment, which results in a range of chemical shifts 20-fold larger than other commonly used nuclei. 2][4][5][6][7] Yet, there remains a large discrepancy between the theoretically predicted and experimentally achieved 129 Xe polarization for continuous-flow production.Even more mysterious, similar operating parameters can yield very different polarization values.Hyperpolarized (HP) 129 Xe is produced using a process called spin-exchange optical pumping (SEOP). 8,9First, the electronic spin of the valence electron of a mediating alkali metal, typically Rb, is polarized using circularly polarized laser light.The electron spin polarization is then transferred to the nucleonic spin of the noble gas isotope, 129 Xe, through the hyperfine interaction during either binary collisions or via the formation of loosely bound Rb-Xe van der Waals molecules. 10,11In this process, the final nuclear spin polarization of 129 Xe depends on the spin-destruction rate, as well as on the Rb polarization and the spin-exchange rate, both of which depend strongly on the Rb vapor density.
For mathematical convenience, current theoretical treatments of continuous-flow SEOP assume that liquid Rb inside the optical cell is in thermal equilibrium with its vapor; even in a closed cell, this condition is only satisfied in the absence of temperature gradients and under specific experimental conditions.In the case of continuousflow SEOP, where gas at room temperature is flown over a pool of hot liquid Rb, the Rb vapor density and distribution inside the optical cell can be quite far from what would be expected under the assumption of thermal equilibrium for a closed cell.Yet, models for homogeneous, temperature-dependent Rb distributions, such as the Smithells and the Killian formula, are commonly used to estimate the Rb vapor density inside the optical cell. 12,13These models do not accurately account for the thermodynamic state of the optical cell, rendering them inapplicable to the continuous-flow SEOP system.
In this paper, by using finite element methods, mechanisms affecting the fluid dynamics that determine the Rb distribution and Xe residence times inside the optical pumping cell, and ultimately the 129 Xe polarization are systematically addressed.In situ integrated atomic absorption spectroscopy measurements are used to measure the Rb density inside the optical cell, thus validating the results of our simulations.

II. METHODS
According to the standard SEOP theory, the polarization of Xe is governed by four parameters: the Rb polarization, P Rb , the spin-exchange rate, γ SE , the spin-destruction rate of Xe, Γ Xe , and the residence time, τ res , In particular, the parameters P Rb and γ SE are known to be strongly affected by the temperature and Rb distribution inside the optical cell.Here, by using finite element analysis, we first determine the temperature and Rb distributions, as well as the Xe residence times for real scale models of optical cells at different volumetric flow rates.We then use these values to create spatial maps of P Rb and γ SE inside the cell.The polarization of each Xe atom is then calculated and iteratively updated every 250 ms, as the gas moves from the inlet to the outlet of the optical cell.Thereby, these simulations account not only for how long the Xe particle spends inside the optical cell, but also for the path taken by the gas.
A. Theory

Effect of Rb density on SEOP parameters
Under optical pumping, the polarization of Rb almost instantaneously enters its steady state given by 10,11 where ε is the fraction of light that is circularly polarized, Γ Rb is the Rb spin-destruction rate, and γ p is the optical pumping rate.Γ Rb comprises of two terms, one from binary collisions and the other from the formation and destruction of van der Waals molecules, 10,11 where and The constant n i represents the density of species i, k i represents the binary cross section of Rb and species i, and b i represents the relative ability of species i to form and break-up van der Waals molecules.The values of k i and b i are reported in Table I.Under most continuous-flow SEOP experimental conditions, Rb-Xe collisions are the dominant Rb spin-destruction mechanism.
The optical pumping rate γ p depends on the local photon flux.The photon flux attenuates down the cell as the light is absorbed by the Rb atoms between the front of the cell and the local point.Therefore, the local optical pumping rate depends on the Rb density TABLE I. Theoretical constants used for polarization calculation.Below are reported the values of the SEOP parameters used in our simulations, with their temperature or field dependence, along with their average values in italics.

Spin-destruction constants of Rb
6:02 Â 10 À21 ( T 298 ) prior to the local point and the overlap of the spectral profile of the laser and the Rb D 1 absorption line.In this work, a 4 cm, 60 W Gaussian beam was simulated to have a wavelength of 794.63 nm and a linewidth of 0.19 nm.The Rb D1 (794.78 nm) had a pressure broadened linewidth of 0.13 nm. 14Assuming the efficiency of circularly polarized light to be unity, the attenuation of the optical pumping rate along the main axis of the optical pumping cell can be described by the following differential equation: [15][16][17] The term β describes the overlap of the Rb absorption profile and the spectral profile of the pump laser, which in our case is 2:61 Â 10 À18 m À2 .The spin-exchange rate, is proportional to the Rb density and comprises of two mechanisms-binary collisions and the formation of van der Waals molecules.The spin-exchange rate due to binary collisions is given by where hσvi, whose value is reported in Table I, is the velocity-averaged Rb-Xe binary spin-exchange cross section.At high gas pressure, binary collisions are the dominant spin-exchange mechanism.As the pressure is lowered, the contribution of spin-exchange due to the formation and destruction of van der Waals molecules increases. 18At pressures at which most continuousflow optical pumping systems operate (4.7 atm in our case), the van der Waals spin-exchange rate is similar to that of binary collisions.Spin-exchange within short-lived van der Waals molecules is more complex than binary collisions as this interaction depends on multiple factors, such as the formation and destruction rates of the van der Waals molecules.Also, as certain atoms are more apt to create these pseudo-Rb-Xe molecules, the spin-exchange rate is also dependent on the densities of the constituent gases and thereby on the gas composition, pressure, and temperature.9][20][21] Here, we just report the result for the Xe, N 2 , and He mixtures used in our experiments, 10,22 where (10) In this expression, κ is the equilibrium constant of the van der Waals molecules, α characterizes the interaction between the alkali valence electron and the nucleus of the noble gas, and Z is the molecular formation rate per unit volume, which sets the magnitude of the spin-transfer rate. 21All of these values are constants for a particular set of atomic species, regardless of the isotopic composition.The constant ζ depends on the specific isotopic composition of Rb and on its polarization.Specifically, in the limit where Rb polarization approaches zero, ζ will assume a value of 0.1791.In the limit where P Rb !100%, as in our case, ζ ¼ 0:0949. 18,20The constants b i describe the adeptness of the third body in creating van der Waals molecules, with a lower b i leading to faster spin-exchange rates for that species. 10,18,20 slight modification of Eq. ( 9), was first found in the original Ruset dissertation in 2005 that was later amended in 2007, 10 This form uses van der Waals specific rates, γ i , for each of the constituent gases, where Equation ( 11) was then adopted by others 2,6,11 and by our group. 17n Ruset, 10 the value of these rates (5230 Hz for Xe, 18 5700 Hz for N 2 , 21,23 and 17 000 Hz for He 24 ) was obtained by using the measurement of Cates of γ M ¼ 2:92 Â 10 4 s À1 , which assumed a Rb vapor density given by the Killian formula, and a value of ζ for the zero Rb polarization limit (ζ ¼ 0:1791).The specific rates of He and N 2 also contain the b i constants 0.31 and 0.92, respectively, while b Xe ¼ 1.Interestingly, b N2 ¼ 0:92 was taken from Rice and Raftery 23 (following a measurement made earlier by Zeng et al. 21) instead of from the measurements made by Cates et al. of b N2 ¼ 0:275. 18An in-depth review of the literature also revealed that an incorrect version of the formula for spin-exchange rate through van der Waals molecules, has been used. 1,16The error led to a misleadingly high spin-exchange rate as individual terms were treated as rates instead of time constants. 22As the accuracy of γ M depends on Rb density, here we will use the γ M value obtained by Shao et al., 25 who precisely measured the Rb density using Faraday Rotation and used the b N2 ¼ 0:275 value measured earlier by Zeng et al. 21Table II compares the difference between the most commonly used parameters, the species-specific rates, the incorrect version that appears in Fink et al. 1,16 and those used for most of our simulations (Table I).As the species-specific rates and our chosen constants (Table I) give rise to spin-exchange rates that differ almost by a factor of 4, both values will be used and compared in the simulations.

Depolarization mechanisms
The means of spin-destruction of 129 Xe comprises two types of interactions: intrinsic and extrinsic. 26The intrinsic interactions consist of the familiar binary collisions and van der Waals molecules, whereas extrinsic interactions consist of the diffusion through magnetic field gradients and wall collisions.A study by Chann et al. showed that because of the weak binding of other Xe gas molecules, the primary source of intrinsic relaxation for densities below 14 amg is from the formation of Xe-Xe van der Waals molecules. 27Further work by Anger led to the formula for intrinsic relaxation, 26,27 where r N2 ¼ 0:51 and r He ¼ 0:25.Wall relaxation due to paramagnetic impurities on the inner surface of the cell is the dominant extrinsic source of spindestruction. 28This means that cell geometry, which will affect the rate of the interaction, is an important factor when determining the wall relaxation rate.Some have suggested a linear dependence on the surface area to volume ratio and a strong suppression of the wall relaxation at high field.These interactions are also known to be temperature dependent, 22,26,27,29,30 which further complicates the matter.A previous study assumed complete Xe depolarization at the cell walls, though they note the atoms will most likely depolarize only after several wall collisions. 16We will speculate the effects of wall collisions on Xe relaxation, but not include them in our simulations due to the large variability noted across the literature.Additionally, we will not consider other spindestruction mechanisms, such as depolarization due to collisions with paramagnetic Rb clusters, that have been proposed to explain the discrepancy between theoretical and experimental 129 Xe polarization values. 6,11fter the Xe leaves the main body of the cell, there are two sources of depolarization that we will consider: depolarization due to dark Rb and solid-state relaxation during cryogenic collection of Xe.As the spin-destruction of Rb is quite high (.10 kHz), it can be assumed that any Rb atom found in the outlet of the optical cell, which is not illuminated by the laser light, is completely depolarized.This dark Rb can still undergo spin-exchange with Xe, thus depolarizing the gas, 17 where P Xe is given by Eq. ( 1).The presence of Rb in the outlet of the optical cell can be inferred from the visible coating of Rb in the u-bend of the outlet of the cell (Fig. 1).In our simulations, we will only include this effect for the time that Xe spends in this region (t out ).For the solid-state Xe relaxation during cryogenic collection, we will consider a collection time t c equal to the one used experimentally and a relaxation time, T 1 , of 87 min. 2When these two sources of depolarization are included, the final Xe polarization will be given by 6,24

B. Finite element method simulations
Simulations of gas dynamics within the optical cell have been done before, for both batch mode and continuous-flow SEOP. 1,16,34ink 16 modeled batch mode SEOP in a cylindrical shaped cell, but assumed a Rb vapor density given by the Smithells formula.This same assumption was made for their simulations of continuousflow SEOP. 1 In this case, while they did model presaturation of the gas mixture from a puddle of Rb placed into a long cylinder, they assumed the gas mixture to be completely saturated with Rb when considering the main body of the cell, again applying an  16 where the different species-specific rates are added together; and (bottom) constants believed to be more applicable for continuous-flow SEOP found in Table I.Note that when the species-specific rates are simply added together, the final spin-exchange rate becomes an order of magnitude larger than the one predicted by the correct formula [Eq.(11)].
γ M (10  empirical formula.Interestingly, they found that a presaturation region of 20-30 cm in length, much longer than the one typically used for % 270 ml cells, would be needed to fully saturate the incoming gas with Rb vapor.Additionally, these simulations have revealed what has been assumed to be turbulence, even for their simplified cell geometries. Here, fluid dynamics simulations were performed on three fullscale models of optical cells (Fig. 2), which are variations of the ones typically used in old models of commercially available polarizers (Polarean 9800, Polarean Inc. Durham, NC, USA).All of the optical cell models had the same main body volume, with an approximate radius of 2.4 cm and a length of 15 cm.They were built in the CAD software SolidWorks (Dassault Systemes SolidWorks Corp, Velizy-Villacoubly, France) and imported into COMSOL Multiphysics (COMSOL, Stockholm, Sweden).Simulations were then performed for a gas mixture consisting of 89% He, 10% N 2 , and 1% Xe at 4.7 atm and for gas flow rates of 0.75, 1.125, and 1.5 SLM.

Boundary conditions
A temperature boundary condition of 20 C was imposed at the inlet of the optical cell (Fig. 1).As polarizers use temperature controllers to keep the cell wall at a given set temperature, another boundary condition of 110 C was enforced on the outer walls of the cell.The surface of the presaturation region was simulated to be at 165 C. The puddle of Rb in the presaturation region was initially modeled to have a surface area of 15.7 mm 2 .To simulate the heat loss from the vaporization of Rb, a modified form of the Hertz-Knudsen equation, was used to quantify the latent heat.Note that the partial pressure of Rb has been removed to simplify computations.Here, H vap is the specific heat of vaporization of Rb, N A is Avogadro's number, α is the sticking coefficient, p sat is the saturation pressure of Rb, m Rb is the molecular mass, k B is the Boltzmann constant, and T is the temperature.The sticking coefficient, which can vary from 0 to 1, was taken to be unity. 1This assumes perfectly efficient evaporation, although with Rb oxidation, this number is expected to decrease over time.As such, the conditions assumed for these simulations assume the best possible scenario for Rb evaporation.The pumping beam contributes significantly to cell heating.In the optical cell laser, heating takes the form of where ν l is the laser frequency and h is Planck's constant. 16This model accounts for the heat capacity of the gas mixture and the Pyrex optical cell and it allows for heat to be dumped into the oven domain as well as carried out of the cell by the flowing gas.

Rubidium distribution
As the optical pumping rate and the laser heating both depend on the Rb distribution, it is important to model how the Rb vapor is being introduced into the optical cell.The evaporation of Rb from a puddle was modeled by the Hertz-Knudsen equation, with the partial pressure of Rb now included, A kinetic model by Ytrehus and Østmo was also tested for the evaporation. 35The diffusion of the Rb vapor throughout the cell was described by Chapman-Enskog theory. 36As the diffusion of Rb is isotropic, it can be treated as though only occurring in the main constituent gas, He.Under this assumption, the diffusion coefficient becomes where P is the pressure and d He and d Rb are the van der Waals diameters of He and Rb.

Preliminary simulations
The first step for our simulations was to calculate the Reynolds number to determine the flow regime.The Reynolds number was determined in three different regions within the optical cell: first at the inlet, then right before the main body of the cell, and finally inside the main body of the cell.The Reynolds numbers were found to be 62, 57, and 66, respectively, for 1.5 SLM, which places the gas dynamics well within the laminar flow regime.Additionally, the Rayleigh number was approximately 2200, which indicates a convection driven flow.As such, the Grashof number, which is analogous to the Reynolds number for convection driven flows, was also calculated.This was found to be 2 Â 10 4 , which again suggests that the flow is laminar as a Grashof number 10 9 indicates a transition to turbulence.To further ensure that a turbulence model was not needed to describe the fluid dynamics, a study was done using the k-ε RANS (Reynolds Averaged Navier-Stokes) model.The parameter that scales how much turbulent kinetic energy is added into the RANS model is known as the turbulent dynamic viscosity, μ T .The maximum value of μ T was of the order of 10 À56 Pa s, which can be compared with the standard viscosity of 10 À5 Pa s.This again indicated the absence of turbulence.As the results of the laminar and k-ε RANS models were identical, further simulations were performed with the laminar flow model, which was computationally more efficient.
Due to the tight coupling of the laser heating to the Rb distribution, both directly and through optical pumping, a fully coupled model that included laser heating could not be solved.Instead, an iterative approach was adopted.First, a coupled model for the fluid dynamics and the thermodynamics was solved using just the first three heat sources discussed: the boundary conditions of the cell walls held at 110 C and presaturation region held at 165 C, and the latent heat of evaporation (study 1).The Rb distribution was then solved for (study 2) using the solution from study 1.A third study (study 3) solved for the optical pumping using the temperature and Rb distributions found in studies 1 and 2. The coordinate dependent laser heating coming from the solution of the optical pumping was then applied to a second model solving for the fluid dynamics and the thermodynamics (study 4), which was followed by another solving for the new Rb distribution (study 5).The solution converged to a stationary state after one iteration of applying the laser heat as described above (one iteration ; studies 1-5).This was determined as the average temperature of the main body of the cell only fluctuated by 0.3% and the average Rb density only fluctuated by 0.7% over three iterations of the described process.Therefore, only one iteration of laser heat was instantiated.

Xe residence times and polarization
Particle tracing simulations were used to calculate the distribution of Xe residence times.For these simulations, 1000 Xe atoms were released from a random distribution right before the main body of the cell.The residence time of each Xe atom was then calculated in COMSOL.Spatial maps of the temperature, Rb density, and Rb polarization were exported to Mathematica (Wolfram Research Inc., Champaign, IL, USA), where the polarization of each Xe atom was iteratively updated every 250 ms.

C. Atomic absorption spectroscopy measurements
In order to experimentally validate the simulations, absorption spectroscopy measurements similar to those that have been previously reported in the literature [37][38][39] were performed to measure the Rb vapor density.These measurements were performed on two different optical cells: a brand-new cell and an older cell.In the latter, the measurements were performed at the peak of the cell's performance (broken-in) as well at the end of its lifetime.For these measurements, a modified oven enclosure for the commercial polarizer was built to allow absorption spectroscopy measurements to be made perpendicular to the beam of the pumping laser.A broadband halogen lamp was used as the light source, similar to Couture et al. 37 A lens, with a 100 mm focal length, was used to collect the light into a fiber optic cable after it passed through the cell.This light was sent into a custom optical spectrometer built in our lab with a resolution of 0.009 nm/pixel.This spectrometer is described in greater detail in Antonacci et al. 17 Spectra were taken immediately after stopping the gas flow.After the absorption spectroscopy measurement, the accumulated Xe was thawed and dispensed into a Tedlar bag (Jensen Inert, Coral Springs, FL, USA) for an absolute polarization measurement made on the Polarean 2281 Polarization Measurement Station.Additional atomic absorption measurements were made on the sealed older cell.

Spectral analysis for Rb density measurements
For these measurements, the intensity of light passing through absorbing atoms can be described by where I o is the characteristic line shape of the lamp, l is the sample length, n is the density of absorbing atoms, and σ is the absorption cross section.By rearranging Eq. ( 20) and integrating over σ, one can extract the density of the species, n, as as by definition, where r o is the classical electron radius, c is the speed of light, and f is the oscillator strength of the transition.The acquired spectra ln(I o =I) were fit with the Lorentzian function, Parameter B accounts for physically irrelevant instrumental gains, A and T represent the depth and asymmetry of the peak, respectively, γ is the peak width, and ν o is the center frequency.The alkali density can then be expressed as where the limits of integration were ten times the fitted width, γ, in order to standardize the analysis across all data sets.This process was repeated for both the D 1 and D 2 peaks, and the results from the two peaks were then averaged.

A. Simulated Rb density and distribution
The Rb density, temperature, and laser heat were found using the iterative method previously described.In addition to the Hertz-Knudsen equation that was used to describe the evaporation of Rb into a flowing gas, a modified kinetic evaporation model by Ytrehus and Østmo was also tested. 35The average Rb density values obtained for the main body of the cell are reported in Table III for various locations of the liquid Rb puddle.Our simulations show a Rb vapor density, originating from Rb puddles of different sizes and locations, that has little to no correlation to the one predicted by empirical saturated vapour pressure curves for alkali metals such as the Killian, the Smithells, or the Nesmeyanov, typically assumed in SEOP theoretical models. 12,13,40As shown in Fig. 3, gas flowing atop of the Rb puddle in the presaturation region is far from becoming saturated.Note that, while for a cell wall temperature of 110 C, the Killian formula predicts a density of 11:0 Â 10 18 m À3 and an even higher density of 2:16 Â 10 20 m À3 for a temperature of 165 C (presaturation region), our simulations predict a Rb density of 7:20 Â 10 18 m À3 for a liquid Rb puddle of 15.7 mm 2 located in the presaturation region.Furthermore, doubling the surface area of the Rb puddle does increase the average Rb density in the main body (Table III).If the puddle of Rb (15.7 mm 2 ) is moved to the main body of the cell (Main), the Rb density drops even further to 1:4 Â 10 17 m À3 .This suggests that Rb puddles in the main body of the cell do not significantly contribute to the overall Rb vapor density during continuous-flow SEOP (Fig. 4).
With use, Rb initially in the presaturation region migrates into the main body of the cell, coating the walls of the main body of the cell.In order to simulate the Rb distribution in a broken-in cell, simulations were run with the inner walls of the cell coated TABLE III.Average Rb densities, main body temperatures, and laser heat deposited at 1.5 SLM for Rb puddles located just in the presaturation region (Presat), just in the main body of the cell (Main), in the presaturation region and in the main body of the cell (Main + Presat), as well as for just a uniform Rb coating of the main body optical cell walls (Coated Walls), and for Rb present both in the presaturation region and on the main optical cell body walls (Coated Walls + Presat).In the first line, the density found using the Østmo kinematic model is noted in parenthesis.The optical cell walls were held at 110 °C and the presaturation coil was held at 165 °C.The laser was assumed to be a 4 cm, 60 W Gaussian beam.with a layer of Rb: first with just the Rb coating on the optical cell walls (Coated Walls), and second with Rb coating of the optical cell walls and a Rb puddle in the presaturation region (Coated Walls+Presat).The results of these simulations can be seen in Fig. 4. Note that when the Rb is in the presaturation region and the cell walls are coated, the saturation vapor pressure does finally almost match the Killian formula, which predicts a density of 11:0 Â 10 18 m À3 , if the cell wall temperature is used.Yet, as seen in Table III, the actual temperature inside the cell is also much hotter than the cell surface and thus, according to the Killian formula, one should expect a much higher Rb density.It is also interesting to note that these simulations show an inhomogeneous distribution of Rb within the optical cell, with a reduced Rb density in the front of the cell, as previously speculated by Appelt et al. 15 Finally, if the cell is closed and contains a puddle of Rb in the main body of the cell, only then does the Rb vapor density reach a value of 10:6 Â 10 18 m À3 .

B. Measured Rb densities
A full table of the experimentally measured Rb densities can be found in Table IV.Rb density measurements were made during different stages of the cell.For the broken-in cell, with Rb present both in the presaturation region and in forming a thin layer of Rb on most of the cell walls, the density measurement indicated a Rb density of approximately 3:1 Â 10 18 m À3 .As the cell aged and polarization started to decline, Rb density dropped by an order of magnitude, down to 2:4 Â 10 17 m À3 .In the new cell (Rb located mainly in the presaturation region, with the main body of the cell having clean walls) measured Rb densities were also low, 6:9 Â 10 17 m À3 , and much smaller than those predicted by the Killian formula for a 110 C cell wall temperature (11:0 Â 10 18 m À3 ).

C. Simulated gas residence time
In continuous-flow SEOP simulations, gas residence time within the main body of the cell is simply calculated using the volumetric flow rate, Q, the total gas density in amagats, [G], and the volume of the cell, V cell by using 6,7 This formula assumes a one-dimensional linear velocity, also known as a plug flow.The results of the velocity field simulations for three different cell geometries, all with the same main body volume, for the fastest and slowest gas flow rate are shown in Fig. 5. Our simulations show a remarkable reduction in residence times with respect to the residence times predicted by Eq. ( 25), which assumes perfectly straight gas pathlines within the cell. 41oreover, our simulations show that there exist different regimes for different flow rates leading to a behavior that is not at all TABLE IV.All density and collected Xe polarization values measured on the broken-in, old, and brand-new optical cells. 129Xe was polarized using a 1/10/89 Xe/N 2 /He gas mixture at the given pressure, temperatures, and flow rates on our Polarean 9800.After an accumulation given by the collect time, Xe was thawed into a Tedlar bag.The polarization was then immediately measured on the our Polarean 2281 Polarization Measurement Station.The closed cell refers to density measurements that were taken on the old cell after it was sealed and removed from the polarizer.Note that in our system, a total laser power of 36 W is measured at the front of the cell.25) suggests.As suggested by the Rayleigh number, the motion of the fluid is dominated by convection, giving rise to two convection cells that are almost symmetric about the primary axis of the optical cell (Fig. 6).These convection cells reduce the effective cell volume, which can be thought of as V eff in Eq. ( 25).These convection cells are also known as dead volume. 41As the flow rate increases, the flow becomes more dominated by the inertial terms, leading to a reduction of the convection rolls, and to an increase in the number of pathlines that move straight through the cell.The distribution of the residence times for all flow rates is shown in Fig. 7.This shows that while some gas particles spend more time in the cell, others will spend less time, as one would expect from mass conservation.A wider distribution indicates that more particles are spending longer times within the cell, which leads to higher overall Xe polarization.A more narrow distribution indicates that more particles are leaving the cell quickly, due to the reduced volume effect, leading to a lower final Xe polarization.This can be seen especially in the flared inlet model as it is shifted left of the other two models, even with similar distribution characteristics (Fig. 7).

D. Simulated Xe polarization values
All parameters needed to calculate the Xe polarization are provided in Table I.We also compared the effect of using the species-specific rates (Table II), previously used in similar work, for the straight inlet model.For these simulations, the Xe polarization was iteratively updated along the path of each Xe atom, using the local temperature and Rb density.As the actual value of Rb density inside the optical cell strongly depends on the size and location of the Rb puddle, to best match our experimental conditions, we used an inflowing Rb density of 3:1 Â 10 18 m À3 , equal to that found experimentally.An example of the Rb and temperature distribution can be found in Fig. 8.Note that while the distribution is slightly inhomogeneous, due to convection bubbles and temperature gradients present in the optical cell, the average density inside the main body is approximately 3:1 Â 10 18 m À3 .
Interestingly, when experimental Rb density values and SEOP parameters reported in Table I were used in our simulations, our simulated Xe polarization values were much lower than what was measured experimentally (Table VI).This was the case even without including in the model other possible depolarization mechanisms that are known to exist, such as wall collisions, depolarization of Xe by dark Rb, or solid-state relaxation during cryogenic extraction.Even more interestingly, when the species-specific rates previously used in the literature for similar calculations 2,6,10,11,17 were used, the simulated Xe polarization dropped even further (Table V).

A. Rb density and distribution
Rb vapor density is the main parameter that affects final Xe polarization.Yet, when modeling Xe polarization values, the Rb vapor density is often deduced from measurements of optical cell wall temperature by using empirical saturated vapour pressure curves for the alkali metal.Our simulations and experimental measurements of Rb density inside the optical pumping cell show that, under experimental conditions typically used for continuous-flow SEOP of Xe, the Rb density inside the main body of the cell is substantially lower than what is predicted by the Killian formula and other empirical formulas.This is in good agreement with previous experiments that found Rb vapor densities to be a factor of 2-3 lower than that predicted by the Killian formula even for closed cells. 25,30While this should not be surprising, as during continuous-flow SEOP, the Rb vapor is far from thermodynamic equilibrium, most theoretical models of continuous-flow SEOP still continue to assume a Rb vapor density close to that of thermodynamic equilibrium. 1,2,6,7,11,17,24,42,43Our simulations also show that the Rb vapor density depends on the size and location of the liquid Rb puddle, and that the condition for fully saturated vapor pressure is often not met for smaller presaturation regions, as the one typically used for small optical cells (cell volume % 270 ml, presaturation region volume % 3:5 ml).Our simulations also show that a few droplets of Rb in the bottom of the main body of the cell do not contribute substantially to the final Rb density under continuous-flow SEOP.Only in the presence of a presaturation region and Rb coating on the cell wall does the Rb density appear to reach an order of magnitude that is comparable to that suggested by the Killian formula.
These simulations also explain why the performance of optical cells improves over time before it starts to decline again.Optical cells are known to have a "breaking-in" period; the Rb density will be much lower than expected in brand-new cells, which are typically ran under hotter conditions than broken-in cells.Polarization remains low until the walls of the main body of the cell are coated with Rb, something that can be easily observed after several polarization runs.This behavior is predicted by our simulations and confirmed experimentally.In very old cells, the presaturation region becomes depleted of Rb, and only the coating on the cell walls remains.This, again, leads to a Rb density much smaller than what one would predict by empirical formulas and what is needed to obtain high polarization values.This unanticipated difficulty in attaining the expected Rb densities may explain why large optical cells that include long Rb presaturation regions 42 historically have performed much better than small cells with short (few cm long) Rb presaturation regions and why the performance of old commercial optical cells, that did not include a Rb presaturation region, was much worse than those of new optical cell models that do include a presaturation region.It is important to note that, to further complicate the issue, with time, Rb oxidation on the cell walls will suppress the amount of Rb that can evaporate, further limiting the Rb vapor density inside the optical cell.
The simulated Rb density appears to reach a density close to that suggested by the Killian formula only in closed cells.This may explain in part why batch mode SEOP systems have outperformed continuous-flow SEOP systems in terms of polarization values achieved. 3These results suggest that, in order to improve optical cell performance, one should use large presaturation regions, while also coating the optical pumping cell walls with Rb.This is typically done for closed cells used in experiments to measure constants relating to SEOP.For instance, Wu et al. 44 baked their closed cell at 80 C for an entire week for the alkali metal and vapor to be as close to thermal equilibrium as possible.It was only then that measured Rb density values agreed with those predicted by the Killian formula.

B. Gas residence times
Xe final polarization depends on Rb density and on the Xe residence time inside the main body of the optical pumping cell.From the simulations of three different full-size optical cell models, it is apparent that residence time depends on the geometry, not just gas flow rate  25) and the corresponding theoretical polarization.The overall shift of the mean residence times from that given by Eq. ( 25) can be attributed to the reduction of the effective volume due to convection rolls.At faster flow rates, particles have enough inertia to quickly leave the cell instead of being trapped in the convection rolls.Particle tracing videos are provided in the supplementary material for a visual representation of gas flowing inside the optical cell.
and volume of the cell.Just by altering the location and style of the inlet, different residence time distributions can be obtained for the same cell volume and flow rate.It is apparent that as the flow rate is lowered, natural convection dominates even more over inertial forces and the distribution of residence times becomes broader (Fig. 7).
Based on the results of these simulations, it is not surprising that better performing polarizers have much longer optical cells. 2,5s extending the main body of the cell is often not an option, as this would require a complete redesign of some polarizers, the solution may lay in simply optimizing the cell geometry.
It is also interesting to note the different behavior of the flow field as a function of flow rate.The standard SEOP theory would suggest that a higher polarization should be achieved at slower flow rates.However, this is not what is experimentally observed in these cells.There appears to be another depolarization mechanism that increases with decreasing flow rates.This has been reported by Norquay et al. 2 as well, who pointed to a study by Walter et al., 45 suggesting that it may be due to laser heating effects at the slower flow rates.By comparing the fastest and slowest flow rates in Fig. 5, it seems as though slower flow rates also increase the number of wall collisions due to the presence of larger convection rolls, possibly leading to a further drop of polarization, as also suggested by Fink et al. 16

C. Xe polarization values
Finally, and more interestingly, when correct Rb densities are used in the simulations, lower Xe polarization values are obtained than those typically predicted theoretically and experimentally observed.This is despite assuming the best experimental conditions and without considering other depolarizing mechanisms that are known to exist, such as wall collisions.One may speculate that the difference between the simulated (% 10% of 60 W) and experimental (% 55% of 36 W) laser absorption could be the origin of such   discrepancy.Note that this difference in laser absorption derives from the distinction between the simulated laser power of 60 W, which assumes no losses, and the true laser power of 36 W, for which an additional reduction in power due to reflections of approximately 15% is expected, 6 and from ignoring Rb wall relaxation.However, as the modeled Rb vapor density was equal to the Rb density measured experimentally and as Rb and Xe wall relaxation was ignored in our simulations, the simulated Rb and Xe polarization values should still be an overestimation of experimental values.Yet, simulated polarization values are actually smaller than those obtained experimentally.This means that there must exist other sources of error in the model.An in-depth analysis of models and parameters previously reported and used in the literature for continuous-flow SEOP of Xe seems to suggest that the issue may be the spin-exchange rate.As current Xe continuous-flow SEOP gas mixtures comprise mostly He, spin-exchange rate is highly determined by b He .The value of b He used in our simulations comes from measurements performed by Ramsey et al. 19 and estimations by Cates et al., 18 where Rb density was deduced from the measurement of the cell wall temperature using empirical formulas.Additionally, these measurements, based on measurements of Xe T 1 , assume a wall relaxation rate that is temperature independent, which recent studies with 3 He suggest that it may not be a correct assumption. 22,30ven more interesting is that, if the species-specific rates are used, the theoretical polarization values are even smaller than what is obtained experimentally.This suggests that, most likely, the Xe spin-exchange rate is much higher than what was previously reported in the literature for lean Xe gas mixtures and under the conditions typically used in commercial continuous-flow SEOP polarizers.Given that some of the key parameters needed to correctly predict final Xe polarization values were measured more than 30 years ago, under different experimental conditions (pressure, temperature, gas mixture, and magnetic field strength) than those used in most continuous-flow SEOP Xe polarizers, it would be advisable to remeasure these parameters, possibly by using methods like those used in Chann for 3 He that are not prone to systematic errors and that do not require accurate knowledge of the Rb vapor density or of the temperature dependence of Xe wall relaxation rate. 30

V. CONCLUSION
In the past, the use of empirical formulas, which assume a thermodynamic equilibrium state for the Rb density inside the optical cell, has artificially inflated theoretical Xe polarization values above experimental ones obtained from continuous-flow SEOP, leading to a search for additional depolarization mechanisms that could explain the discrepancy between theoretical and experimental Xe polarization values. 6ere, more realistic fluid dynamic simulations of continuousflow SEOP dynamics, along with atomic absorption spectroscopy measurements, indicate that the Rb vapor density inside the optical pumping cell during continuous-flow SEOP is far below what one would expect under the assumption of thermodynamic equilibrium.When these more realistic Rb densities are used in the standard SEOP model, theoretical Xe polarization values drop below those obtained experimentally, even when optimal conditions are assumed.This discrepancy suggests that some of the parameters previously used in the standard SEOP model, in particular, the spin-exchange rate, may not be accurate and will need to be remeasured under the same conditions used for continuous-flow SEOP of similar Xe gas mixtures.

SUPPLEMENTARY MATERIAL
See the supplementary material for particle tracing videos of the straight inlet model at 0.75 SLM and 1.5 SLM.

FIG. 1 .
FIG. 1. CAD model of the optical pumping cell with a straight, bottom inlet.Rb is typically placed in the presaturation region.With use, a visible Rb coating is found on the walls of the main body of the cell and in the u-bend portion of the cell outlet, indicated in the figure.

FIG. 3 .FIG. 4 .
FIG. 3. Rb vapor distribution from a Rb puddle in the presaturation region shown as a slice through the cell for a flow rate of 1.5 SLM: (top) main body of the cell and (bottom) the presaturation region.

FIG. 5 .
FIG. 5. Velocity fields for the slowest (0.75 SLM) and fastest (1.5 SLM) flow rates.Note the larger number of convection rolls present at 0.75 SLM vs 1.5 SLM.

FIG. 6 .
FIG. 6. Streamlines at 1.5 SLM in the straight inlet cell.Two large convection rolls symmetric about the main axis, which decrease the effective volume of the cell, are evident.

FIG. 7 .
FIG. 7. Xe residence time distribution plot for the different optical pumping cell geometries and flow rates simulated as well as the corresponding Xe polarization distribution.A wider distribution, shifted to the right, consistently produces better polarization.Dotted vertical lines indicate the Xe residency times calculated by using Eq.(25) and the corresponding theoretical polarization.The overall shift of the mean residence times from that given by Eq. (25) can be attributed to the reduction of the effective volume due to convection rolls.At faster flow rates, particles have enough inertia to quickly leave the cell instead of being trapped in the convection rolls.Particle tracing videos are provided in the supplementary material for a visual representation of gas flowing inside the optical cell.

FIG. 8 .
FIG. 8. Rb density in m À3 (left) and temperature in C (right) for straight inlet cell at 1.5 SLM with an inflowing Rb density of 3:1 Â 10 18 m À3 , wall and presaturation coil temperatures of 110 C and 165 C, respectively, pressure of 4.7 atm, and a 60 W laser.The gas mixture was 1/10/89 of Xe/N 2 /He.

TABLE II .
Comparison of the density normalized spin-exchange rate due to the formation and breakup of van der Waals molecules for (top) the most common species-specific rates; (middle) the incorrect version found inFink et al.,

TABLE V .
Comparison of polarization after the main body using the constants listed in TableIand the previously used species-specific rates for the straight inlet model.

TABLE VI .
Simulated polarization values for each flow rate, for all models, right after the main body of the cell [Eq.(1)], accounting for dark Rb in the u-bend of the outlet [Eq.(15)], and finally accounting for solid-state relaxation (T 1 = 87 min) during a 25 min collection of Xe in the frozen state [Eq.(15)].Note that the depolarization due to dark Rb is negligible in most cases.The flared inlet, which is the model closest to that currently on our polarizer, is also compared to the experimental polarization, which is reported as the average measurement for the broken-in cell (boldface values).