Analytical solution of light diffusion and its potential application for light simulation in DUNE Vyacheslav Galymov IPN Lyon DUNE DP-PD Consortium Meeting 02.11.2017
Some challenges S2 Electron • Light simulation for dual-phase has to include drift time • Generation of S2 in addition to S1 • Light conversion on the cathode plane if used S1 • The challenging aspect is how to populate PMTs with a photons produced along particle tracks TPB/ITO coated cathode • The solution so far to produce a light map (or (Not an option?) light library in larsoft) which defines visibility of a given detector voxel wrt to the photon detectors • Note: time spread due to RS is not applied to PMT array photon arrival times in larsoft • Size of the map can quickly become a Since we are not interested in tracing challenge due to large detector volume paths of each photon, but rather the • Simulation of light visibility from each voxel, end result, is it possible to find an although to be done once, also becomes a effective theoretical description? CPU intensive task 2
Photon transport in diffusion media • Actually there has been a big interest in this question due to its medical applications to evaluate light propagation in tissues (e.g., oxygen meters) • Also in nuclear physics: neutron transport O 2 meter(image: Wikipedia) Basically find effective solution for particle propagation in scattering medium using diffusion theory 3
Diffusion equations • Generally described by Fokker-Plank (FP) PDE: 𝜖𝑢 𝑞(𝑦, 𝑢) = 𝐸 𝜖 2 𝜖 𝜖 𝜖𝑦 2 𝑞 𝑦, 𝑢 − 𝑤 𝑒 𝜖𝑦 𝑞(𝑦, 𝑢) Where is D is constant diffusion coefficient and 𝑤 𝑒 is constant drift velocity • For 𝑤 𝑒 = 0 FP PDE reduces to differential equation describing Brownian motion (Wiener process): 𝜖𝑢 𝑞(𝑦, 𝑢) = 𝐸 𝜖 2 𝜖 𝜖𝑦 2 𝑞 𝑦, 𝑢 This is the equation one needs to solve for photon diffusion subject to appropriate boundary conditions 4
Boundary conditions Photon are absorbed on the cathode absorption condition for this plane For other sides of the TPC, the simplest assumption is that photons exiting TPC do not contribute in any significant way absorption boundary would also be appropriate But could also consider a quasi-reflective boundary at some point Absorption boundary condition: p=0 𝑞(𝑦, 𝑢) 𝑇 = 0 Reflective boundary condition: p=0 𝑞(𝑦, 𝑢) 𝑇 = 𝑑𝑝𝑜𝑡𝑢 𝜖 𝜖𝑦 𝑞(𝑦, 𝑢) 𝑇 = 0 5
Diffusion from a point source In unbound medium solution for diffusion equation for point source at 𝑠 0 , 𝑢 0 is given by Green’s function: 𝒔 − 𝒔 0 2 1 𝐻 𝒔, 𝑢; 𝒔 0 , 𝑢 0 = 3/2 exp − 4𝜌𝐸𝑑 𝑢 − 𝑢 0 4𝐸𝑑 𝑢 − 𝑢 0 Where c is the velocity of light in the medium. For LAr c = 21.7 cm/ns 1 𝐸 = 𝜈 𝐵 - absorption coefficient [1/units of L] 3(𝜈 𝐵 + (1 − )𝜈 𝑇 ) 𝜈 𝑇 - scattering coefficient [1/units of L] – average scattering cosine For 𝜈 𝑇 = 1 55 and 𝜈 𝐵 ~0 Isotropic scattering = 0 • • Including Ar form factors introduces some 𝐸 = 18.8 cm anisotropy for Rayleigh scattering = 0.025 Or cm2/ns if one multiply by velocity to get more familiar units 6
Unbound solution Time profile for source 3m away from detector Note the extending tail is due to infinite boundaries due to scattering photons will keep arriving … 7
Single absorption boundary Time profile for source 3m away from detector Inf boundary 𝐻 ∞ Image S True S 𝐻 𝐶 𝑦 0 𝑏 2𝑏 − 𝑦 0 Solution for 𝑦 > 𝑏 is simply a difference between two unbound Green’s functions The tail is reduced due to photons for true source 𝑦 0 at and its mirror absorbed at the boundary image at 2𝑏 − 𝑦 0 𝑞 𝑦, 𝑦 0 , 𝑢 = 𝐻 𝑦, 𝑦 0 , 𝑢 − 𝐻(𝑦, 2𝑏 − 𝑦 0 , 𝑢) 8
Source between two absorbing planes Source b/w two absorption boundaries at -a and a Image S - Image S + True S - 𝑏 𝑦 0 𝑏 2𝑏 − 𝑦 0 Could use image source method as well, but need to also absorb image sources at further boundary: in the sketch that would be S − (−2a + 𝑦 0 ) at boundary 𝑏 would need an image source at 4𝑏 + 𝑦 0 and so on Just like an image of a mirror reflection in a mirror or a screen capture of a screen capture on a video call Of course each contribution becomes smaller and smaller correction truncates the infinite series 9
Source reflection Reflection operations: • Negative boundary at -a: -2a – x • Positive boundary at +a: 2a – x First few terms in the series Image source Add/Subtract Img Source 1 Img Source 2 −𝑦 ′ − 2𝑏 −𝑦 ′ + 2𝑏 1 - 𝑦 ′ − 4𝑏 𝑦 ′ + 4𝑏 2 + −𝑦 ′ − 6𝑏 −𝑦 ′ + 6𝑏 3 - … … … … Subtract terms with n/2 = odd, add terms with n/2 = even 10
Full solution 1D 𝜖𝑢 𝑞(𝑦, 𝑢) = 𝐸 𝜖 2 𝜖 Diffusion PDE: 𝜖𝑦 2 𝑞 𝑦, 𝑢 with absorption at x ± 𝑏 +∞ exp − 𝑦 − 𝑦 ′ + 4𝑜𝑏 2 − exp − 𝑦 + 𝑦 ′ + 4𝑜 − 2 𝑏 2 𝑞(𝑦, 𝑢) ∝ 4𝐸𝑢 4𝐸𝑢 𝑜=−∞ 11
Solution for point source in 3D 𝜖𝑦 2 𝑞 + 𝜖 2 𝜖 2 𝜖𝑧 2 𝑞 + 𝜖 2 𝜖 𝜖𝑢 𝑞 = 𝐸 𝜖𝑨 2 𝑞 With absorbing boundaries at 𝑦 𝑐 = ±𝑥 , 𝑧 𝑐 = ±𝑚 , 𝑨 𝑐 = ±ℎ , Take: 𝑞 = 𝑌 𝑦, 𝑢 × 𝑍 𝑧, 𝑢 × 𝑎(𝑨, 𝑢) 3D PDE reduces to 1D PDE for each component 2 𝑌 𝜖 𝑢 𝑌 = 𝜖 𝑦 Since 1D has been solved, we have simply to 2 𝑍 𝜖 𝑢 𝑍 = 𝜖 𝑧 take a product of 1D solutions 2 𝑎 𝜖 𝑢 𝑎 = 𝜖 𝑨 12
Full solution in 3D 1 𝑞 𝒔, 𝑢; 𝒔 0 , 𝑢 0 = 3/2 × 𝑇 𝑦 × 𝑇 𝑧 × 𝑇 𝑨 4𝜌𝐸 𝑢 − 𝑢 0 +∞ exp − 𝑦 − 𝑦 0 + 4𝑜𝑥 2 − exp − 𝑦 + 𝑦 0 + 4𝑜 − 2 𝑥 2 𝑇 𝑦 = 4𝐸(𝑢 − 𝑢 0 ) 4𝐸(𝑢 − 𝑢 0 ) 𝑜=−∞ +∞ exp − 𝑧 − 𝑧 0 + 4𝑜𝑚 2 − exp − 𝑧 + 𝑧 0 + 4𝑜 − 2 𝑚 2 𝑇 𝑧 = 4𝐸(𝑢 − 𝑢 0 ) 4𝐸(𝑢 − 𝑢 0 ) 𝑜=−∞ +∞ exp − 𝑨 − 𝑨 0 + 4𝑜ℎ 2 − exp − 𝑨 + 𝑨 0 + 4𝑜 − 2 ℎ 2 𝑇 𝑨 = 4𝐸(𝑢 − 𝑢 0 ) 4𝐸(𝑢 − 𝑢 0 ) 𝑜=−∞ This gives us photon concentration density in any point at any given time 13
Source at (0,0,0) in a 6x6x6 box t = 100 ns t = 10 ns Infinite solution Bounded solution t = 1000 ns Particles have diffused to the walls where they were absorbed 14
Photon flux across the surface What is of interest to us is the so-called time of first passage The time photon hit a given surface The overall integral of this distribution would give us an acceptance probability for this point Note that by construction 𝑞 𝒔, 𝑢 𝑇 = 0 Fick’s law of diffusion relates flux to the concentration density: 𝐾 𝒔, 𝑢; 𝒔 0 , 𝑢 0 = −𝐸𝛼𝑞(𝒔, 𝑢; 𝒔 0 , 𝑢 0 ) The change in particle density crossing the surface per unit time: 𝜖 𝑢 𝑄 Ω 𝑢; 𝑠 0 , 𝑢 0 = 𝒆𝑩 ∙ 𝐸𝛼𝑞 Ω 15
Photon flux PDF at a bounding surface 3D PDF in the volume: 1 𝑞 𝒔, 𝑢; 𝒔 0 , 𝑢 0 = 3/2 × 𝑇 𝑦 × 𝑇 𝑧 × 𝑇 𝑨 4𝜌𝐸 𝑢 − 𝑢 0 And the Cartesian components of the flux vector are 𝑘 + 𝑇 𝑦 𝑇 𝑧 𝜖 𝑨 𝑇 𝑨 𝐾~𝑇 𝑧 𝑇 𝑨 𝜖 𝑦 𝑇 𝑦 𝑗 + 𝑇 𝑦 𝑇 𝑨 𝜖 𝑧 𝑇 𝑧 𝑙 Since we are working with a cubical geometry the unit 𝑘 , ± normal to each face would simply be ± 𝑗 , ± 𝑙 So depending on the face the integrand 𝒆𝑩 ∙ 𝑲 reduces to one of a the appropriate J term 16
Photon flux PDF at a bounding surface Consider we are interested at surface z = -300 (e.g., cathode plane in 6x6x6) 1 𝑔 𝑦, 𝑧, 𝑢; 𝑦 0 , 𝑧 0 , 𝑨 0 , 𝑢 0 = 3/2 × 𝑇 𝑦 × 𝑇 𝑧 × 𝜖 𝑨 𝑇 𝑨 4𝜌𝐸 𝑢 − 𝑢 0 𝑨=−300 Derivative wrt z Independent of z now But still depend on of z0 evaluated at z = -300 Since we have a sum of Gaussians of the form 𝐻~ exp[−𝑡 𝑦 − 𝑦 0 2 ] 𝜖 𝑦 𝐻 = −2𝑡 𝑦 − 𝑦 0 𝐻 17
Integration Spatial integrals can be done quickly For the acceptances calculation need to integrate Gaussians in the expansion series of the type 𝑐 𝑒𝑦 exp[−𝑡 𝑦 − 𝑦 0 2 ] ~ erf … 𝑏 Interpolate error function table computed in advance fast and independent of integration range, since only need two end-points 𝑦 Tabulate 0.5 erf 𝜏 2 up to N𝜏(= 1) = 𝑂 • 𝑦+𝐸 𝑦 Integral for any interval [x+D, x] 𝑇 − 𝑇 • 𝜏 𝑦 𝜏 𝑦 𝑏 𝑐 • Integral for an interval [-a, b] 𝑇 𝜏 𝑦 + 𝑇 𝜏 𝑦 18
Acceptance calculation: basic sanity check This gives the acceptance per 𝑒𝑢 𝒆𝑩 ∙ 𝐸𝛼𝑞 detector face Ω For a cubical boundary and the source at the center the answer is simply : 1/6 ≈ 1.666667 Calculation gives exactly that! More detailed comparison can be done against MC simulation of photon transport 19
Recommend
More recommend
