A Novel, Finite-Element Model for Spin-Exchange Optical Pumping Using an Open-Source Code
AA Novel, Finite-Element Model for Spin-Exchange Optical Pumping Using anOpen-Source Code
G.M. Schrank (Dated: September 14, 2020)A new model is presented for spin-exchange optical pumping using an open-source code,ElmerFEM-CSC. The model builds on previous models by adding the effects of alkali-vapor hetero-geneity in optical pumping cells and by modeling the effects of hyperpolarized-gas wall-relaxationusing a diffusion model. The code supports full, three-dimensional solutions to optical-pumpingmodels, and solves for (1) laser absorption, (2) alkali vapor concentration, (3) fluid flow parameters,(4) thermal effects due to the pumping laser, and (5) noble gas polarization. The source code forthe model is available for researchers to utilize and modify.
I. INTRODUCTION
Spin-exchange optical pumping (SEOP) is a techniquewhereby the ensemble nuclear spin-angular momentumof certain noble gasses can be increased to of order 10%.The technique is currently most notably used clinicallyand pre-clinically in lung imaging using MRI [1], but ithas also been used in the NMR characterization of porousmedia [2] and protein dynamics [3].The physics of SEOP are described comprehensivelyin other places [4, 5]. Briefly, the technique involves twosteps: (1) optical pumping of an alkali-metal vapor and(2) spin-exchange from the alkali-metal vapor to a noble-gas nuclei. In the first step, optical pumping, a beamof circularly-polarized light is directed onto a transpar-ent cell containing a macroscopic amount of alkali metal.The cell is heated, usually to between 100-200 ◦ C, in or-der to vaporize some amount of alkali metal. The laserinteracts with the metal vapor to create close to 100%spin polarization of the alkali vapor.In the second step, spin-exchange, the alkali vaportransfers spin-angular momentum to the noble gas nu-clei. A gas mixture is introduced into the cell containinga noble gas and some other inert gasses, and through col-lisional interactions, the metal vapor transfers its spin-polarization to the noble gas. The alkali metal vaporatoms become depolarized in this interaction, but be-cause of optical pumping, the alkali atoms are quicklyrepolarized.One popular method of SEOP involves the hyperpo-larization of xenon-129 (
Xe) using rubidium (Rb) va-por. Hyperpolarized
Xe (HP
Xe) gas is typicallyproduced in a continuous manner using a flow-throughpolarizer [6–8]. A flow-through polarizer operates byflowing a
Xe gas mixture through an optical pump-ing cell containing Rb vapor. The
Xe spin-exchangeinteraction occurs on a short enough timescale that con-siderable
Xe polarizations can be achieved during theshort transit through the optical pumping cell.Currently, there are no SEOP models that attempt toaccount for the full three-dimensional flow-dynamics ofoptical pumping cell geometries. Computational modelsfor SEOP can be broken into two groups: finite differencemodels and finite element models (FEM). The finite dif- ference models appear to have been the first in use, butdetailed descriptions of the models do not appear exten-sively in the literature. Models of this type approximateflow through the cell by either a one-dimensional, plug-flow model (such as used in Ref. [9]) or two-dimensional,laminar-flow model(such as Ref. [10]).The first computational model of spin-exchange opticalpumping to be extensively described in the literature wasdone so by Ref. [11]. It originally described a model ap-proximating an optical pumping cell geometry with onlya half-cylinder. The model included the effects of laserheating, fluid dynamics, and heat transfer. The originalmodel only investigated a static optical pumping cell inwhich the gas was not flowing. This was later extendedin Ref. [12] to analyze two geometries with gas flow:a half cylinder and a simplified geometry of the opticalpumping cell used by Ref. [7].Recently, Ref. [13] has described an FEM fluid modelof an optical pumping cell. However, it does not appearthe model included SEOP, and it only included fluid flow,heat transfer, and diffusion of Rb metal vapor.Here, we present an open-source FEM code thatattempts to model the dynamics of SEOP using thefull three-dimensional geometry of a common optical-pumping-cell design. The model incorporates fluid flow,Rb vapor diffusion, thermal transfer, laser absorption,and
Xe polarization.
II. MODEL DESCRIPTION
The FEM model presented here utilizes the finite-element method to numerically solve five differentialequations. The author used an open-source code,ElmerFEM-CSC, to implement the model. The modelsimulates (1) fluid flow through the cell, (2) diffusion andtransport of the Rb vapor, (3) heat transfer through thegas stream, (4) absorption of the laser by the Rb andsubsequent polarization, and (5) spin-transfer betweenthe Rb vapor and the
Xe. Of the five modeled quanti-ties, ElmerFEM-CSC contained modules for solving thefirst three. Descriptions of these three modules can befound in Ref. [14]. The last two equations, for laser ab-sorption and
Xe polarization, were obtained by modi- a r X i v : . [ c ond - m a t . o t h e r] S e p fying existing ElmerFEM-CSC modules to accommodatethe expressions that describe the dynamics of those pro-cesses.That this is an open-source model is a notable dif-ference from the other models described in section I.The other FEM models were produced using commercialcodes that are not easily accessible to many researchers.Other finite-difference codes have not been formally pub-lished in an open-access environment.The specifics of the expressions used to describe laserabsorption and Xe polarization very closely followedRef. [11]. Major differences between the current modeland other models will be highlighted in the following sec-tions.
A. Laser Absorption Model
Laser absorption is one of the key features of an SEOPmodel. The expression for modeling laser absorption bythe Rb-metal vapor is [15]: ∂ψ ( ν, z ) ∂z = − n Rb ( z ) σ s ( ν ) Γ SD ( z ) γ opt ( z ) + Γ SD ( z ) ψ ( ν, z ) (1)where ψ is the photon flux density, z is the azimuthal spa-cial coordinate, ν is the frequency of the light, n Rb is thenumber-density of the Rb, σ s is the cross-section for ab-sorption by unpolarized Rb, Γ SD is the spin-destructionrate of the Rb, and γ opt = (cid:82) ∞ ψσ s ∂ν is the optical pump-ing rate. The expression is easily solved using finite-difference methods, and it was used in many previousmodels. However, the presence of the integral expressionis challenging for standard finite-element methods.Instead, the FEM model presented here uses themethod described by Ref. [11]: ∂γ opt ∂z = − βγ opt n Rb (cid:18) − γ opt γ opt + Γ SD (cid:19) (2)where β = 2 (cid:112) π ln(2) r e f λ l w (cid:48) (cid:16) i (cid:112) ln (2)( r + is ) (cid:17) δλ . (3)Here, r e is the classical radius of the electron, f is theoscillator strength of the Rb D-1 line, λ l is the laser wave-length, and δλ l is the laser-line width. w (cid:48) ( Z ) is the realpart of the Faddeeva function; with s denoting the ratioof the laser-to-Rb-line-frequency difference and the laserline width, and r denoting the ratio of the Rb-absorption-line width and the laser-line width.Equation (2) solves explicitly for γ opt with the assump-tion that the beam’s spectral profile is Gaussian through-out the optical pumping cell. This is notably differentfrom finite-difference models, in which the shape of thespectral profile changes during passage through the opti-cal pumping cell. The laser absorption solver used the ElmerFEM-CSCAdvection-Reaction module as its template. The equa-tion solved by the Advection-Reaction Module is [14]: ∂c∂t + (cid:126)v · (cid:126) ∇ c + Γ c = S (4)The following modifications were made to eq. (4).First, (cid:126)v can be constrained to be a unit vector, ˆ n ,pointing in the direction of the laser beam propagation.The spatial derivative of eq. (2) can be rewritten asˆ n · (cid:126) ∇ γ opt = ∂γ opt ∂z .Second, Γ from eq. (4) can be set equal to:Γ = − βn Rb (cid:18) − γ opt γ opt + Γ SD (cid:19) . (5)The non-linear portion of the equation is solved itera-tively by the Picard method [14].The source term, S , is set to 0. B. Diffusion of the Alkali Metal
In many previous models, the Rb-metal vapor distri-bution was assumed to be uniform. In this FEM model, adiffusion model of the Rb-metal vapor is implemented us-ing the diffusion module supplied with ElmerFEM-CSC.Simulated geometries include temperature-dependentsources and sinks, which use the Hertz-Knudsen equa-tion for the boundary condition of the source/sink [12]: j Rb = α Rb p sat ( T ) − p √ πM Rb k B T (6)where p sat ( T ) is the saturation partial pressure for agiven temperature, p is the instantaneous partial pres-sure, M Rb is the molecular mass of the Rb, k B is Boltz-mann’s constant, T is the absolute temperature, and α Rb is the evaporation coefficient of Rb. The saturation par-tial pressure p sat ( T ) for Rb is calculated by the Killianequation [16]. The evaporation rate j Rb can be positive(sources) or negative (sinks) depending on the local val-ues of p sat ( T ) and p .Although α Rb has not been measured for Rb, Ref. [12]notes that the ideal value of α Rb = 1 is expected. Forall the simulations presented here, this value of α Rb wasused. For the current model, the saturation and instan-taneous partial pressures were converted to the absolutemass concentration used as default in ElmerFEM-CSCdiffusion module [14]. C. Wall-Relaxation of HP Xe In most previous simulations, wall-relaxation is mod-eled as a constant term in the HP
Xe spin-relaxationterm. In this simulation, the diffusion-based model ofHP
Xe wall-relaxation is a refinement of the expressionpresented in Ref. [11]. In that model, the wall boundary-conditions were modeled as completely depolarizing HP
Xe spin-polarization at the walls. The authors offeredan alternative model with the depolarization set to 1%rather than 100%. However, they noted that the lack ofexperimental data hindered more precise estimates.The present model attempts to refine this approxima-tion and connect the boundary condition at the walls tothe wall-relaxation time. The wall-relaxation time (atroom temperature) can easily be measured for a givencell by filling the cell with HP
Xe and monitoring theamplitude of the HP
Xe NMR as a function of time.It is known that this relaxation time is typically 10s ofminutes [9].The present model uses the solution to the diffusionequation for surface-evaporation in a sphere. Like withthe evaporation of the Rb, we wish to have an expressionthat relates the flux of polarization to the surface: j = − D Xe ∂P∂r = α RW P, (7)where α RW is the so-called mass transfer coefficient.The solution to the diffusion equation for a spherical FIG. 1. A plot of the decay curve (eq. (9)). The modeledsphere had a radius of R =0.03 m and a diffusion coefficient of D Xe = 1 ∗ − s . The mass transfer coefficient of α RW =3 . ∗ − was calculated using eq. (11) for a wall-relaxationtime of τ wall =300 sec. The overall decay curve (solid blue)overlaps with the first-order term (dashed red). The second-order term does not contribute more than 1 part per thousandto the overall decay curve, and it decays to less than 10 − afterless than 100 seconds. The higher-order terms contribute lessand decay even more quickly. This wall-relaxation time iseven shorter than discussed in the text, and illustrates therobustness of the assumption at reasonable wall-relaxationtimes. Compare the overall decay curve with the decay curvefrom Figure 2 . FIG. 2. The visualization of an FEM diffusion model using thesame parameters as those in used in Figure 1. The left graphicdepicts the concentration of the species (normalized to one)at τ wall = 300 sec. The graph on the right is the averagedconcentration as a function of time. Note that the curve isqualitatively identical to the theoretical curve in Figure 1. geometry with this boundary condition is [17]: P ( r ) P i = 2 LRr ∞ (cid:88) n =1 exp( − D Xe β n t/R ) β n + L ( L −
1) sin( β n r/R )sin β n (8)where the β n ’s are the roots of β n cot ( β n ) + L − L = Rα RW D Xe , P is the polarization after time t , P i is theinitial polarization, and R is the radius of the sphere.The total polarization, P tot , decay curve is given by [17]: P tot P i = ∞ (cid:88) n =1 L exp( − D Xe β n t/R ) β n [ β n + L ( L − . (9)The wall-relaxation time, τ wall , is usually found by fit-ting the decaying HP Xe NMR amplitude to an expo-nential function: P tot P i = exp (cid:18) − tτ wall (cid:19) . (10)In the limit of β n >> β ; ∀ n (cid:54) = 1, the larger β n s can beignored, and only the β term will significantly contributeat long timescales. In this case, we can compare the time-dependent portion of just the first term of eq. (9) to eq.(10), and we find that τ wall = R D Xe β .It turns out that for typical values found in SEOPsystems (i.e. R > D Xe ≈ . cm s and, τ wall > β ≈
1, in which case
L < .
5, the ratio β β (cid:39) β n s are much larger. For wall-relaxationtimes longer than 10 min., the ratio between the β n s and β is even greater.Using this approximation, we can solve for α RW andfind that: α RW = D Xe R − (cid:115) R D Xe τ wall cot (cid:115) R D Xe τ wall . (11)Equation (11) is a poorly behaved equation and needsto be applied with care. In particular, for a given τ wall ,the expression has asymptotes at R = πn √ D Xe τ wall ; n = 1 , , , ... . The spacing of the asymptotes is propor-tional to multiples of the diffusion length. Physically,this indicates that the HP Xe cannot undergo wall re-laxation faster than it takes for a polarized
Xe atom atthe center of the sphere to travel to the wall. Therefore,eq. (11) is only valid for τ wall > R π D Xe .As an example of the validity of the assumption thatthe higher-order terms in eq. (8) and (9) do not signifi-cantly contribute, Figure 1 shows the theoretical contri-butions to the decrease in concentration of the polarizedspecies as a function of time for a sphere. The first-orderterm almost precisely overlays the overall decay. Thehigher-order terms quickly decay and do not significantlycontribute to the decay curve.The expression from eq. (11) was used to model thewall-relaxation of HP Xe that diffused to the wall. AnFEM model using this expression was checked computa-tionally in a simple spherical model to give the correcttransient relaxation time and decay curve (see Figure 2).
D. Thermal Transfer through Cell Walls
The present model attempts to capture the nuancesof thermal transfer through the optical pumping cellwalls. The model assumes the optical pumping cell isin a forced-air oven where the temperature of the air isheld constant. Heat transfer in the walls of the cell issolved by imposing the boundary condition: − k ∂T∂n = α T ( T − T ext ) (12)where k is the heat conductivity of the gas mixture inthe cell, T is the temperature at the boundary, T ext isthe temperature at which the external flowing air is held,and α T is the heat transfer coefficient. The heat transfercoefficient is calculated by combining the affects of theoptical pumping wall (usually glass; thermal conductivity k wall = 1 . Wm · K ) and that of forced-air convection (heattransfer coefficient α air ≈ Wm · K ). These two termscan be combined into a single heat transfer coefficientboundary condition by using the expression [18]:1 α T = t wall k wall + 1 α air (13)where t wall is the thickness of the wall. The expressionassumes that the contact area of the different boundariesare approximately equal, which is true when t wall is smallcompared with the other linear dimensions of the optical-pumping cell. E. Other Considerations
Other quantities, such as viscosity and heat capacity,were calculated using standard expressions (see Table I). These quantities were used as inputs for the standardmodules in ElmerFEM-CSC. For more information, seeappendix A.
III. VERIFICATION
Although the model contains the relevant SEOPphysics, it is necessary to verify that the model is provid-ing an adequate computational solution to those expres-sions. In order to confirm that the FEM model conformedto accepted results, a comparison study was run with themodel described in Ref. [9], the Freeman model.The Freeman model differs from the model describedhere in several important ways. The Freeman model usesa finite-difference method rather than a finite-elementmethod. The Freeman model is a one-dimensional sim-ulation. The fluid velocity, temperature, and the Rb-density distribution are approximated as uniform in theFreeman model. The laser spectral profile used in theFreeman model is not constrained to be Gaussian; ratherthe spectral profile is discretized and calculated at eachspacial node in the model. Finally, in the Freeman model,the
Xe polarization is calculated using the average Rbpolarization rather than Rb polarization at each spacialnode.It was possible to reproduce aspects of this model inthe infrastructure of the FEM model described here. Thegeometry of the FEM model was drawn as a cylinder.The Navier-Stokes fluid model, heat transfer model, andRb-diffusion model were disabled. They were replacedwith uniform axial flow, a uniform temperature, and ahomogeneous Rb number-density. In addition, for thiscomparison, the wall-relaxation described in section II Cwas replaced by an additional term in the HP
Xe re-laxation expression. The outlet
Xe polarization of thisFEM model was compared with the predicted
Xe po-larization from the Freeman model.The results of this comparison are shown in Figure 3.The two models predicted very close to the same polar-ization at low temperature for all flow velocities. Thediscrepancy between the models increases with both thetemperature and the flow velocity. The Freeman modelnearly always predicts a lower polarization. The largestabsolute and relative discrepancies are 9 percentage-points (16%) and 19% (8 percentage-points), respec-tively.A couple of explanations for this discrepancy were con-sidered and rejected. One explanation for this discrep-ancy is that the calculated laser absorption differs in thetwo models. The Freeman model does not enforce the as-sumption of a Gaussian spectral profile of the laser, whilethe FEM model described here does. The assumption of aGaussian spectral profile might influence the calculatedRb polarization, and thus the final
Xe polarization.However, when a comparison of the 0.2 SLM flow-rateseries was conducted, it was both found that the averagecalculated Rb polarization differed by no more than 0.2
FIG. 3. The predicted
Xe polarization from the Freemanmodel and the FEM model as a function of temperature andflow velocity using substantially similar assumptions regard-ing flow, temperature distribution, and Rb number-densitydistribution. Each curve represents a different flow velocity(0.2, 0.5, 1.0, and 1.5 SLM) from a different model. The dot-ted lines are the predictions of the Freeman model, and thesolid lines are the predictions of the FEM model. The Free-man model almost always predicts a lower polarization thanthe FEM model. This is thought to be because the FEMmodel includes a combination of a diffusion model for the
Xe polarization solver and a heterogeneous Rb polariza-tion spatial-distribution. percentage-points (0.2%) and that the FEM model usu-ally predicted a lower
Rb polarization. Thus, it seemsunlikely that the Gaussian spectral profile, which is en-forced by the FEM model, would result in a higher pre-diction for the
Xe polarization.Another possible cause for the discrepancy between thetwo models was the Rb polarization spacial-distributionthat existed in FEM model. Although the Freemanmodel does calculate a spacial-distribution for the Rbpolarization, the final
Xe polarization calculation usesa spatially-averaged value. In order to test if this was thecause of the discrepancy, the flow direction in the FEMmodel was reversed so that it was directed parallel to thedirection of laser propagation instead of anti-parallel.The Rb polarization in the FEM model is highestwhere the laser first enters the cell geometry. There-fore, reversing the direction of the flow in the modelchanges the spatial-distribution of the Rb polarizationfor the simulated
Xe. In the anti-parallel flow con-figuration, the gas will first encounter Rb with a lower-than-average Rb polarization and progress to a region ofhigher-than-average Rb polarization. In the parallel flowconfiguration, this order is reversed.Although the reversed-flow FEM model did result ina lower prediction of the
Xe polarization, this affectalone was not sufficient to account for the discrepancybetween the Freeman model and the FEM model. Thepredicted
Xe polarization at 150 ◦ C from the FEMmodel decreased by only ∼ FIG. 4. A visualization of the geometry chosen on which totest the FEM code. The geometry is a 100-cc standard opticalpumping cell. The entire geometry (optical pumping region,inlet, and outlet) was drawn using Onshape CAD software,and a mesh was created from this geometry using Salome9.2.1. The Rb source was modelled as a thin, cylindrical filmthat encircled the middle of the optical pumping region ofthe cell (black section in figure). An additional Rb sink wasprescribed in the entire outlet tube for the geometry (notpictured). the FEM model incorporates diffusion of the
Xe po-larization while the Freeman model does not. This differ-ence, coupled with a Rb polarization spatial-distributiondiscussed above, will cause the FEM model to predicthigher
Xe polarizations than the Freeman model bysignificant amounts. The solution to the one-dimensional diffusion-reaction equation , assuming a linear Rb polar-ization as a function of axial position, causes the finalpredicted
Xe polarization to fall off linearly as a func-tion of gas velocity. Solving the advection-reaction equa-tion , which is used in the Freeman model, with a Rb po-larization independent of axial position causes the finalpredicted
Xe polarization to fall off exponentially withgas velocity. The details of this derivation are providedin appendix B.This final combined effect from both the diffusion-reaction equation and the Rb-spatial distribution is likelythe cause of the discrepancy between the predictions ofthe Freeman model and the FEM model.
IV. TEST ON A STANDARD SEOP GEOMETRY
After verification, limited tests were conducted on astandard, SEOP-cell geometry to assure that all mod-ules would work in concert. A 100-cc-cell geometry withRb in the optical-pumping body was chosen as the firstgeometry to test (see Figure 4). The test was conductedwith the following boundary conditions: • Laser Power: 75 W • Laser Spectral Width: 0.3 nm • Laser Beam Radius: 1 mm smaller than opticalpumping region radius
FIG. 5. A visualization of the 100-cc cell simulation. The im-age on the left shows the predicted flow-lines for the gas as itmoves through the cell. The flow lines are color-coded by thepredicted temperature of the gas at each point. The slice ofthe cell is color-coded by the predicted Rb number-density inthe cell. The graph on the right shows three of the simulatedquantities averaged over the optical-pumping region of thecell: (1) Rb polarization, (2) absorbed laser power, and (3)temperature. The Rb polarization and absorbed laser powerare both normalized to one. The temperature is normalizedto the fraction above the set-point temperature, in this case,110 ◦ C. • Wall Temperature: 110 ◦ C • Flow Rate: 1.5 SLM • Wall-Relaxation Time: 56 min. • Fraction
Xe: 1% • Fraction N : 10% • Fraction He: 89% • Cell Pressure: 73 psig (assumed to be at sea level)Initially, a steady-state solution for the model was at-tempted. However, when this failed to converge, themodel was solved as a transient problem with 0.1 sec.per time-step. A total of 611 time-steps were simulated.However, due to a file-corruption error, the last 6 time-steps were discarded.Due to possible instability with the module, the
Xepolarization module was not used during the tran-sient simulation. Instead, after the transient simulationreached steady-state in (1) temperature, (2) Rb polar-ization, and (3) laser power absorbed, the last transientsolution was used as the initial conditions for a steady-state calculation with the HP
Xe polarization moduleenabled. Steady-state was defined as less than a 0.05%change in all of the three other metrics.The transient simulation shows convergence to asteady-state solution after the 611 time-steps (Figure 5).The predicted
Xe polarization is comparable to theobserved
Xe polarization, made in Ref. [9], from a100-cc cell similar to the geometry simulated in the FEMmodel. In the Ref. [9] study, a 100-cc cell was observedto absorb 30% of the incident light and produce HP
Xepolarized to ∼
15% at 90 ◦ C. The FEM model predictedthat the simulated 100-cc cell would absorb 7.9% of theincident light and produce HP
Xe polarized to ∼ FIG. 6. A visualization of the steady state solution to theHP
Xe polarization. The visualization on the left showsthe polarization in the optical pumping section of the cell.The maximum polarization is ∼ at 110 ◦ C. The HP
Xe polarization for the FEM simu-lation was taken as the average polarization across a slice5 cm from the edge of the outlet of the model’s geometry.The discrepancy between the observed light absorp-tion of the 100-cc optical pumping cell and the simulatedcell may be from a couple sources. First, as previouslymentioned, the FEM model enforces a Gaussian spec-tral profile throughout the optical pumping cell. Thisassumption may affect the calculated total absorption ofthe light. The affects of the imposed Gaussian distribu-tion are not entirely understood, however, as stated insection III, the Gaussian spectral distribution was foundto not significantly alter the average Rb polarization inthe FEM model when compared with a model that doesnot enforce a Gaussian spectral distribution.Second, the Rb-source distribution in the model maynot accurately describe the distribution in the opticalpumping cell tested in Ref. [9]. The authors do notdescribe the Rb-metal distribution in their 100-cc cell(e.g. a Rb droplet, a thin layer of Rb on the side of thecell, etc.). The distribution and location of the Rb-metalin the cell may affect the details of the dynamics in theoptical-pumping region.The discrepancy between the observed HP
Xe po-larization may be similarly explained by the Rb-sourcedistribution in the model. Because the Rb-source distri-bution may not reflect the conditions in the actual cellthat was tested, the Rb number-density calculated in themodel may be lower than the number-density that wasrealized in the actual cell. This would result in a lowerspin-exchange rate than what was realized in the actualcell, and thus, a lower predicted HP
Xe polarization.
V. SUMMARY AND FUTURE WORK
A new open-source, FEM-based SEOP model has beencoded, tested against an existing model, and tested witha three-dimensional geometry. The model provides pre-dictions of HP
Xe polarizations that are comparable toexisting, accepted models when the conditions are lim-ited to the scope of those existing models. The three-dimensional model predicts HP
Xe polarizations thatare comparable to existing observations. The new FEMmodel provides the ability to visualize important SEOPphenomena such as laser heating, Rb-vapor distribution,and gas flow. The model can compute solutions for com-plicated geometries including current designs for opticalpumping cells.The most immediate improvement that can be madeis to add the ability to multi-thread computations. Thelack of multi-threading limits the speed at which solu-tions can be computed, which in turn, limits transientstudies to simulate only the first several hundred secondsafter initialization. The three-dimensional simulation de-scribed in section IV required ∼ ∼
60 sec. of the cell after initial-ization. Important, long-term behavior cannot be thor-oughly investigated. Multi-threading capabilities withappropriate increases in other computational resourcescould potentially decrease the amount of time requiredto calculate a single time-step and open up the possibil-ity of exploring the long-term behavior of these systems.Better approximations to spin-exchange parametersand other model parameters may increase the fidelity ofthe model’s predictions. However, it is suspected that amajor source of error in the model is likely due to thechoice of Rb source distribution. Because all the param-eters are correct to an order of magnitude, it is unlikelysubstantial gains will be realized by pursuing more accu-rate approximations to these parameters unless it can beshown that the current estimates of the parameters areincorrect by large margins.Further simulations on various optical pumping cell ge-ometries will be reported in a follow-up report.
ACKNOWLEDGMENTS
The author gratefully acknowledges J. Cook for use-ful conversations regarding the coding of the simulations.Also, the author is grateful for insightful discussions withand feedback on initial drafts from B.J. Anderson, T.Barthel, B. Driehuys, and B. Saam. Finally, much of thecomputation time on Amazon’s AWS Cloud ComputingService was supported by backers of the project from Ex-periments.com. The author gratefully acknowledges thefinancial support of all the backers of the project.
Appendix A: Model Details
Section II was meant to highlight the important dif-ferences between the FEM model and previous models.However, readers that are interested in using the modelmay be interested in some more details of that model.The most recent version of the model can be found at https://github.com/drschrank/elmerfem and is freefor use under the GNU license requirement found in thesource code.Although the model can be used in a steady-state for-mulation of the differential equations, transient simula-tions are frequently useful because initial simulations onsome optical pumping geometries fail to converge whensolving the steady-state equations. During testing, thetransient time-step for the simulation was 0.1 sec/step.Two different linear solvers with different convergencelimits were used for the modules. The generalized mini-mal residual method (GMRES) and biconjugate gradientstabilized method (BiCGSTAB) were both used for theNavier-Stokes equation. The BiCGSTAB was used ex-clusively for all the other modules. Non-linear and linearconvergence limits were all less than 10 − and were usu-ally 10 − .As described in section II, the model consists of fivemodules which are calculated in sequence. The firstmodule which is calculated is the solution to the Navier-Stokes equation. It was provided with ElmerFEM-CSCand was used without modification. ElmerFEM-CSCprovides two applicable compressiblity models for gasflow: the “Perfect Gas” compressibility model, whichsolves the full Navier-Stokes equation, and the “In-compressible” compressibility model, which solves theNavier-Stokes equation with the assumption that (cid:126) ∇ · j =0. The “Incompressible”compressiblity model requiresless computational resources, and it is appropriate whenthe Mach number of the fluid flow is low, which is the casefor all reasonable SEOP models. When using the “In-compressible” model, the density can be calculated as afunction of temperature at each time-step. The Navier-Stokes calculation is dependent on the solution to theheat equation, which is discussed later.The second module which is calculated is Rb diffusion.This module uses the Advection-Diffusion module. Thismodule was also provided with ElmerFEM-CSC and wasused without modification. This module’s solution is onlydependent on the resulting gas flow calculated by theNavier-Stokes equation.The third module is the laser absorption module. Thismodule was based on the Advection-Reaction (eq. (4))module, which was provided by ElmerFEM-CSC. Themodule was modified to use the Picard method to ap-proximate the non-linear solution by a series of iterativesteps (see Ref. [14] section on linearization of the Navier-Stokes equation for an example). The specific equationsolved using this method is described in section II A. Thesolution to this module is only dependent on the solutionto the previous module, the Rb diffusion module.The time-dependent portion of eq. (4) is ignored byElmerFEM-CSC solvers when the model is used to solvesteady-state problems. However, for the transient simu-lations, the time-dependent term had to be included. Tohandle this, the time-dependent term was multiplied bya coefficient that is much smaller than the characteristictime-steps used in the model. This modification effec- TABLE I. The table lists all of the expressions that are used to calculate various parameters in the model. The first columnlists the term in the model. The second column lists in which model the term is specifically used. The third column liststhe expression used. Note that because the model is coded in SI units, the actual implementation of some of the expressionmay have been multiplied by a constant to make all the units consistent. The final column lists the reference from which eachexpression is derived. In all cases, the notation used by the reference has been kept. For the meaning of the particular notationin the expressions, the reader should refer to the particular reference. *Note: A previous version of this pre-print incorrectlylisted the expressions used for specific heat ratio and viscosity.
Term Model Expressions Reference
Gas MixtureDensity Navier-StokesHeat ρ = PTR R = C p γ − γ [14]Gas MixtureViscosity Navier-Stokes µ = µ (cid:48) (cid:0) TT (cid:48) (cid:1) T (cid:48) + ST + S µ tot = (cid:80) i x i ( µ i ) Heat Capacityat Constant Pressure Navier-StokesHeatRb Diffusion C p,mix = (cid:80) i x i C p,i [19]Specific HeatRatio Navier-StokeHeatRb Diffusion γ = (cid:80) i x i γ i Alkali EvaporationRate Rb Diffusion j Rb = α p sat − p √ πM Rb k B T [12]Diffusion Constant Rb Diffusion Xe Polarization D , = . × − T / (cid:113) M + M pσ , Ω [18]Optical PumpingRate Laser Absorption ∂γ p ∂z = − βγ p n Rb (cid:16) − γ p γ p +Γ SD (cid:17) [11]Alkali SpinDestruction Rate Laser Absorption Γ SD = Γ Rb + Γ Xe + Γ N + Γ He + Γ V W Γ Rb = κ Rb [ Rb ]Γ N = 170 (cid:16) T − ◦ C . ◦ C (cid:17) [ N ]Γ He = 24 . (cid:16) T − ◦ C . ◦ C (cid:17) [ He ]Γ Xe = κ Xe [ Xe ]Γ V W = f Xe +1 . f N +3 . f He [10][20][21][22]Xenon Spin-Exchange Rate Xe Polarization γ SE = κ SE [ Rb ] [11]Xenon Spin-Relaxation Rate Xe Polarization Γ SR = Γ B + Γ V W Γ b = κ Xe [ Xe ]Γ vdW = Γ XevdW r [ B ][ Xe ] [23]Wall-Relaxation Xe Polarization α = D Xe R (cid:16) − (cid:113) R D Xe τ cot (cid:104)(cid:113) R D Xe τ (cid:105)(cid:17) [17]Thermal Conductivity Heat k = 1 . × − √ T/Mσ Ω k mix = (cid:80) α x α k α (cid:80) β x β φ β [18]Heat TransferCoefficient Heat U = h + (cid:80) nj =1 x j − x j − k j − ,j − h n [18]Laser Heating Heat Q = hν l n Rb γ p Γ SD γ p +Γ SD [11] tively causes each time-step to be a steady-state solutionof eq. (2).The fourth module is the Xe polarization module.This module is based on the Advection-Diffusion modulewhich was provided by ElmerFEM-CSC. The method bywhich the equation was solved was not changed. Theonly changes made were the assignments to the variousconstant parameters of the Advection-Diffusion equation,and the solver was forced to always use the “absolutemass” setting because the equation in this form is easilyadapted to eq. (B1).The solution to this module is only dependent on thelaser absorption module, and the solution does not effectany of the other modules. This allows for the possibilityof running simulations with only the other four modulesactive and then using a final, steady-state solution of theother four modules to calculate the solution to the
Xepolarization.The final module is the heat equation module. Thismodule was provided by ElmerFEM-CSC without anymodification. The solution is dependent on both the laserabsorption solution and the Navier-Stokes solution.Parameters for the various modules are listed in tableI. The values of particular constants used in the equationsin the table can be found in the references listed in thetable.
Appendix B: Details of the Derivation of theDifference Between the Freeman Model and theFEM Model
In section III, it is stated without proof that the as-sumptions of the Freeman model give rise to an expo-nential dependence on flow rate when a uniform Rb po-larization is assumed, while the assumptions of the FEMmodel give rise to a linear dependence on flow rate whena Rb polarization gradient is assumed. In this appendix,the derivation of that result will be given.The FEM model uses the advection-diffusion equationto model
Xe polarization. The form of that equationin one dimension is: D Xe ∂ P Xe ( x ) ∂x + v ∂P Xe ( x ) ∂x +( σ + Γ) P Xe ( x ) = σP Rb ( x ) . (B1)The Rb polarization is assumed to have a linear gradi- ent given by: P Rb ( x ) = P L − P L x + P (B2)where P L is the highest polarization of the Rb at thepoint x = L , and P is the lowest Rb polarization at thepoint x = 0. From Ref. [4], it is clear in the limit as L → ∞ that P Xe → σσ +Γ P L . The solution to eq. (B1) isgiven by: P Xe ( x ) = K e x v + √ v − DXe ( σ +Γ)2 DXe + K e x v − √ v − DXe ( σ +Γ)2 DXe − vσ ( P L − P ) L ( σ + Γ) + σ ( P L − P ) xL ( σ + Γ) + σP ( σ + Γ) . (B3)For v >> D Xe ( σ + Γ), this can be simplified to: P Xe ( x ) = K e − vxDXe + K + σσ + Γ (cid:18) P Rb ( x ) − v ( P L − P L ( σ + Γ) (cid:19) . (B4)If vLD Xe >>
1, then P Xe ( L ) = K + σσ + Γ (cid:18) P L − v ( P L − P L ( σ + Γ) (cid:19) . (B5)For L → ∞ , we get: P Xe ( L → ∞ ) → K + σσ + Γ P L = σσ + Γ P L , (B6)which implies: K = 0 . (B7)Therefore, solving the one-dimensional advection-diffusion equation for the polarization of Xe at x = L assuming a linear Rb polarization gradient gives: P Xe ( L ) = σσ + Γ (cid:18) P L − v ( P L − P L ( σ + Γ) (cid:19) , (B8)which linearly decreases as the flow rate, v , increases.The Freeman model uses eq. (B1) with the diffusionterm absent and an average Rb polarization, ¯ P ave = P L + P . The solution to this equation can be triviallyshown to be: P Xe ( L ) = ¯ P ave σσ + Γ (cid:16) − e − ( σ +Γ) Lv (cid:17) , (B9)which exponentially decreases as a function of v . [1] A.-M. Oros and N. J. Shah, Physics in Medicine andBiology , R105 (2004).[2] V. V. Terskikh, I. L. Moudrakovski, S. R. Breeze, S. Lang, C. I. Ratcliffe, J. A. Ripmeester, and A. Sa-yari, Langmuir , 5653 (2002).[3] L. Schrder, Physica Medica , 3 (2013). [4] T. G. Walker, Journal of Physics: Conference Series ,012001 (2011).[5] S. Appelt, A. B.-A. Baranga, C. J. Erickson, M. V. Ro-malis, A. R. Young, and W. Happer, Phys. Rev. A ,1412 (1998).[6] B. Driehuys, G. D. Cates, E. Miron, K. Sauer, D. K.Walter, and W. Happer, Applied Physics Letters ,1668 (1996).[7] I. C. Ruset, S. Ketel, and F. W. Hersman, Phys. Rev.Lett. , 053002 (2006).[8] G. Schrank, Z. Ma, A. Schoeck, and B. Saam, Phys. Rev.A , 063424 (2009).[9] M. S. Freeman, K. Emami, and B. Driehuys, Phys. Rev.A , 023406 (2014).[10] I. C. Ruset, Hyperpolarized
Xe Production and Ap-plications , Ph.D. thesis, University of New Hampshire(2005).[11] A. Fink, D. Baumer, and E. Brunner, Phys. Rev. A ,053411 (2005).[12] A. Fink and E. Brunner, Appl. Phys. B , 65 (2007).[13] A. Burant, Characterizing Hyperpolarized
Xe Depolar-ization Mechanisms During Continuous-Flow Spin Ex-change Optical Pumping and as a Source of Image Con-trast , Ph.D. thesis, University of North Carolina, Chapel Hill (2018).[14] P. R˚aback, M. Malinen, J. Ruokolainen, A. Pursula, andT. Zwinger,
Elmer Models Manual (CSC IT Center forScience, 2015).[15] M. E. Wagshul and T. E. Chupp, Phys. Rev. A , 4447(1989).[16] T. J. Killian, Phys. Rev. , 578 (1926).[17] J. Crank, The Mathematics of Diffusion , 2nd ed. (OxfordUniversity Press, London, 1975).[18] R. B. R. B. Bird, W. E. Stewart, and E. N. Lightfoot,
Transport phenomena (J. Wiley, 2007) p. 905.[19] “NIST Chemistry WebBook NIST Standard ReferenceDatabase Number 69,” (2018).[20] I. A. Nelson,
Physics of Practical Spin-Exchange OpticalPumping , Ph.D. thesis, University of Wisconsin-Madison(2001).[21] D. K. Walter, W. M. Griffith, and W. Happer, Phys.Rev. Lett. , 093004 (2002).[22] A. Ben-Amar Baranga, S. Appelt, M. V. Romalis, C. J.Erickson, A. R. Young, G. D. Cates, and W. Happer,Phys. Rev. Lett. , 2801 (1998).[23] B. Chann, I. A. Nelson, L. W. Anderson, B. Driehuys,and T. G. Walker, Phys. Rev. Lett.88