Predicting Transport Effects of Scintillation Light Signals in Large-Scale Liquid Argon Detectors

Liquid argon is being employed as a detector medium in neutrino physics and Dark Matter searches. While argon scintillation light has been the primary observable in low-energy dark matter experiments, neutrino detectors have predominantly used it for trigger purposes. A recent push to expand the applications of scintillation light in Liquid Argon Time Projection Chamber neutrino detectors has necessitated the development of advanced methods of simulating liquid argon scintillation light. The presently available methods tend to be prohibitively slow or imprecise due to the combination of detector size and the amount of energy deposited by beam neutrino interactions. In this work we present a semi-analytical model to predict the quantity of vacuum ultra-violet argon scintillation light observed by a light detector, based only on the relative positions between the event and light detector with a precision better than $10\%$. We also provide a method to predict the distribution of arrival times of these photons accounting for the effects of Rayleigh scattering, absorption and reflections. Additionally, we present an equivalent model to predict the number of photons and their arrival times in the case of a wavelength-shifting, highly-reflective layer being present on the detector cathode. Our proposed method can be used to simulate light propagation in large-scale liquid argon detectors such as DUNE or SBND.


I. INTRODUCTION
Several experiments setting out to search for the elusive dark matter particles [1,2] or perform precision measurements of neutrino parameters [3][4][5] have chosen liquid argon (LAr) as their detector medium. Liquid argon is relatively dense (1.41 g/cm 3 ) and is chemically inert allowing ionisation charge to be drifted for distances of several metres. These properties combined with its relatively low price make liquid argon an excellent choice for large scale detectors [6], enabling the construction of modules as large as several kTons [5]. In addition to the ionisation charge used by Liquid Argon Time Projection Chamber (LArTPC) neutrino detectors, liquid argon is an excellent scintillator emitting on the order of 40000 photons per MeV of deposited energy. This scintillation light has been employed by experiments searching for dark matter for energy reconstruction and background rejection [1,7]. In neutrino detectors however, liquid argon scintillation light has not been as thoroughly exploited to date. The MicroBooNE [3] and ICARUS [8] experiments have mostly used scintillation light as a means of triggering and rejecting cosmic events. It has been recently proposed [9][10][11] that LArTPC neutrino detectors with enhanced light collection capabilities could employ scintillation light for improved time, calorimetry and position resolution.
More sophisticated applications of scintillation light in LArTPC neutrino detectors will require precise simulation of the light to develop new algorithms and val- * Corresponding author: dgarciag@ugr.es idate their performance. Such simulation quickly becomes computationally challenging when the detector size reaches tens of tons (or even kTons as in the case of DUNE) combined with the events of interest depositing hundreds of MeV of energy, as expected for accelerator neutrinos. Simulating each photon emitted by the interacting events individually becomes prohibitively slow. A solution applied in the LArSoft software package [12] used by most LArTPC detectors is to implement a lookup library method [13]. In this method the computationally challenging work is performed once, by running a large stand-alone job that generates a vast number of isotropic photons from pre-defined cubes (voxels) inside of the detector active volume. For all voxel-photon-detector pairs the probability of light being detected is saved in a dedicated file. This file is then accessed by standard simulation jobs that use it to estimate the number of detected photons for a given energy deposition without having to simulate each photon separately. This approach works reasonably well in estimating the amount of light for detectors of tens of tons such as MicroBooNE or SBND, albeit introducing a certain granularity into the simulation driven by the size of the voxels used. However, for extremely large detectors this approach results in very large lookup files of several GB even when considering relatively large voxel sizes. Additionally, to develop applications of scintillation light involving time, a good understanding of the effects of light propagation in the medium is needed in addition to the understanding of the intrinsic scintillation light components [14] and the effects of wavelength shifting, e.g. by tetraphenylbutadiene (TPB) [15]. At the large distances present in currently built and proposed LArTPCs, Rayleigh scattering becomes an important factor. This scattering can result in a smearing of the photon arrival times leading to non-trivial effects in their arrival time distributions. The lookup-library method is not designed to predict these effects and, given the non-trivial distributions, a bruteforce approach could necessitate storing a function definition (or a histogram) for each voxel-photon-detector pair. Incorporating this could greatly increase the size of the required lookup-library.
We propose a method to numerically predict the number of photons observed in a particular photon-detector based only on the size and position of the energy depositions in the liquid argon volume. This approach allows fast simulation of scintillation light in large scale liquid argon detectors with a precision better than 10%, even for cases of large numbers of photons originating from high energy particle interactions. In addition, we have developed an analogous method to predict the number of detected photons arriving from from a wavelengthshifter coated, highly-reflective, detector cathode. Passive light detection elements of this type, in the form of reflective TPB-coated foils, are planned to be installed in the SBND experiment [16] and have been proposed as an option for the DUNE detectors [17]. Finally, we also present a method to generate the distribution of photon arrival times that accounts for effects of Rayleigh scattering in liquid argon, and reflections and absorption on the detector walls, using only the relative positions of the energy deposition and the photon detectors.
This article is structured as follows: we first present the basics of liquid argon scintillation light emission and detection, that are relevant to the simulation method we propose. We then describe the framework used to perform the studies described in this work. In Section IV we describe the model to predict the number of photons arriving at a detector given only the position and scale of the energy deposition for both direct light as well as light reflected off the cathode of the TPC. In Section V we describe the model to predict the distribution of arrival times of the photons for both of the above cases. We then test the performance of these models compared to a standalone Geant4 simulation, and to predictions obtained using lookup libraries. Finally, we show an example of the application of this model to a realistic event.

A. Production of Scintillation Light in LAr
The liquid argon scintillation light originates from the de-excitation of argon dimer states formed when an argon atom, excited or ionized by an interaction of a charged particle, attaches a neutral argon atom. This results in a relatively narrow emission wavelength with a peak of 128 nm [14], in the vacuum ultra-violet (VUV) range. This mechanism is very prolific resulting in over 40000 photons being emitted per MeV of deposited energy in the absence of an electric field [18]. Applying an elec-tric field reduces the recombination of argon dimers, decreasing the scintillation yield. At an electric field of 500 V/cm, a typical value used in LArTPCs, the scintillation yield decreases to about 20000 photons/MeV [19].
The time distribution of the emitted photons is composed of two exponential decaying functions with lifetimes of ∼6 ns and ∼1.5 µs, corresponding to the two possible dimer molecular states: the singlet and the triplet [20].
The amount of light emitted can decrease through quenching effects (Q) or recombination (R). Quenching either caused by the ionisation density [21,22] or because of contaminants [23,24] that can absorb the energy of the de-excited dimer without emitting light. Recombination depends on the value of the electric field E in the argon as well as on the ionisation density.
The ratio between the amount of light emitted by each of the two components depends on the ionisation density created by the interacting particle, and is therefore used as a method of particle identification, e.g. in Dark Matter experiments [25]. The decay times, particularly of the slow component, can also be affected by contaminants present in the argon such as nitrogen [23] and oxygen [24] through a quenching process where some dimers transfer their energy to the contaminants instead of de-exciting. The time composition can also be altered by doping with other noble gases, for example xenon [26].

B. Transport of Scintillation Light in LAr
The scintillation light emitted by the argon dimers is able to travel long distances in the liquid argon. Its mean free path is primarily affected by Rayleigh scattering and absorption on contaminants. Rayleigh scattering does not change the number of travelling photons but deflects them on their path. This can be either detrimental or beneficial for the probability of light arriving at photondetectors, depending on the position in the detector and distance travelled. Light undergoing scattering and still arriving at a photon-detector will have travelled a longer path than light impinging on the detector without any interactions. This leads to a non-trivial distribution of arrival times of the photons that could be interpreted as an apparent lengthening particularly of the fast component of the scintillation light. The value of the Rayleigh scattering length (λ RS ) is currently under intense study with different measurements and theoretical predictions reporting values from around 50 cm [27][28][29] all the way up to 110 cm [30]. In this paper we use a value of 100 cm which is close to the most recent reported value [31], and inside of the range of the other expected values.
Argon scintillation photons can also be absorbed (Q abs ) by contaminants with a high cross-section for VUV photons such as nitrogen [32] and other elements that have been observed in commercial argon [28]. The total absorption can be modelled as an exponential suppression of the number of photons as a function of the  [31]. A line at 128 nm is drawn to guide the eye to the scintillation emission wavelength in argon.
distance travelled with the absorption length as a parameter. Due to the high refractive index of liquid argon at VUV wavelengths [33] the group velocity of photons emitted at 128 nm is about two times slower than that of light at visible wavelengths, see Fig. 1.

C. Detection of Scintillation Light in LAr
Detecting argon VUV scintillation light most often requires photon-detectors (PD) able to operate at, or close to, liquid argon temperatures. Cryogenic Photomultiplier tubes have been the default technology used in liquid argon detectors for some time [34,35] and have reported Quantum Efficiencies (QE) of up to 30%. Recently, the idea of using Silicon Photomulipltiers (SiPM) has been gaining ground due to their low power consumption, small size, excellent noise performance at liquid argon temperatures and high QE up to 40% [36]. SiPMs can be used as-is or enhanced using light-collectors such as light-guide bars [37] or a light-trap, such as the ARA-PUCA or X-ARAPUCA devices [38]. Another crucial challenge of detecting LAr scintillation light is the registration of VUV light before it is absorbed by the materials, e.g. glass or plastic, used to shield the sensitive area of the PD. The solution most commonly employed is to coat the PDs with a wavelength-shifting (WLS) com-pound, which absorbs the VUV light and emits light in the visible spectrum easily detectable by the PDs. The direction of the re-emitted light is random, so WLScoated PDs suffer a ∼50% decrease in efficiency due to light emitted away from the active surface.
The travel time, t t , of the scintillation light in large liquid argon detectors ranges from a few to several tens of nanoseconds. Wavelength-shifting compounds tend to impact the photon arrival time, because the emission of the visible light is not instantaneous and has an intrinsic decay time, t W LS , which could delay the detection of the photons. The electronics and data acquisition chains in LArTPC detectors are usually designed with a resolution in a similar range, from 1-2 ns to 16 ns in most recent liquid argon neutrino experiments [3,5,39].
In large LArTPCs used for neutrino experiments the PDs are usually placed behind planes of sense wires [3,8,10]. These and other components of the detector can introduce a further decrease in the number of detected photons due to shadowing effects (Q trans ).

D. Passive elements of detection, i.e. reflective foils
The majority of materials used to construct LArTPC detectors absorb VUV photons causing them to be lost. A method to recover them used primarily in direct dark matter search detectors is to cover the walls of the detector with a highly reflective surface, e.g. ESR foils [1] or PTFE [7], coated with a wavelength-shifting material. These surfaces become passive elements of the light detection system and enhance the amount of light detected by the PDs. It should be noted that the light arriving at the PDs is now already shifted to visible wavelengths where the efficiency of the WLS-coated PDs could be different. Additionally, the group velocity and Rayleigh scattering length of photons at visible wavelengths in LAr are significantly different, see Fig. 1, which will have an impact in any transport effects.
Recently, neutrino experiments have begun exploring the method of installing WLS-coated reflector foils on the detector cathode. Examples include the LArIAT experiment (fieldcage and cathode) [10] and SBND [16]. Implementing this solution has also been proposed for the DUNE detectors [17].

E. Generic Model for Predicting Behavior of Scintillation Light Photons
In general, the number of photons detected by a given PD from an energy deposit ∆E at position (d, θ) can be calculated using the formula: where the scintillation yield S γ (E ) = R(E )/W ph is the number of photons emitted per unit of deposited energy at an electric field E . This is defined in terms of the work function, W ph = 19.5 eV [18], or average energy needed to create a photon at E = 0, and the recombination factor, R(E ), that accounts for the reduction of the scintillation yield due to the presence of an electric field. The position (d, θ) is defined in terms of the distance, d, between the energy deposit and the PD, and the offset angle, θ, between the energy deposit and the normal to the PD surface. Q is the quenching at emission, Q det is the PD efficiency, Q abs (d) is the loss due to absorption effects and depends on d, and Q trans (θ) represents the loss of transmission due to shadowing effects and depends on θ. All of the Q x parameters have values in the range from 0 to 1. P (d, θ) is the geometric coverage of the PD and T (d, θ) signifies other transport effects, both of which depend on distance and angle. Most of the parameters in Eq. 1 are independent of each other and during simulation can be applied at the stage that the light is generated, with the exception of P and T that are tied together and more complicated in their application. The focus of this paper is a method to estimate these two quantities. Similarly, the time at which a photon is detected by a particular PD can be obtained by summing together the independent components resulting from the different processes that photons undergo: where t E is the emission time determined from a distribution combining the dimer decay times with any quenching of the time components. t t (d, θ) is the transport time resulting from the different paths the photons can take to arrive at the PD. t W LS is the time resulting from the intrinsic relaxation time of the wavelength-shifting compound. Finally, t det is the time due to the PD and electronics response. These components can be applied separately and independently. In this paper we also propose a model to calculate t t (d, θ).

III. SIMULATION FRAMEWORK
To develop, test and validate the model of scintillation light transport described in this work, we compared it against a Geant4 [40] simulation embedded in the LAr-Soft software framework [12]. Geant4 is capable of simulating liquid argon scintillation light emission, transport and boundary properties. We use these results as a baseline (i.e. our true information). In a Geant4 simulation, the user needs to provide the optical properties of the active medium, liquid argon, and all of the surrounding materials with which the photons can interact. The optical properties of LAr that we implemented are summarized in Table I. Additionally, to account for potential contaminants in the detector we apply an absorption length of λ abs = 20 m corresponding to 3 ppm of nitrogen [32].
We define our simulated detector geometries using the GDML markup language [43] to be realistic models of real-life LArTPC detectors. Figure 2 shows a schematic  [42] refractive index spectrum 1.32 [31] Rayleigh scattering length spectrum 100 cm [31] TABLE I. Liquid argon properties used in the simulation.

Field-Cage
Active volume PD-plane representation of our geometry. The main volume of liquid argon is delimited by metal walls of a cryostat 1 . The field-cage surrounding the active volume is modelled as an array of metal strips on the top, bottom, upstream and downstream walls and spaced to provide ∼30% optical transparency (a typical value in real detectors). The PDs are uniformly distributed in the open plane of the field-cage (left), referred to as the photon-detector plane (PD-plane). We simulate the PD sensitive windows as flat disks or rectangles facing the active volume. The cathode plane (right) is modelled as an opaque volume covered by polymetric reflector foils coated in a WLS. The WLS used in our simulations is TPB, for which the absorption/emission spectra are taken from [44]. The reflectivities of the materials to VUV photons (pure scintillation emission) and visible photons (coming from the re-emission of the VUV photons absorbed by the TPB on the reflector foils) are listed in Table II. In the simulations used in this work, the reflections on optical boundaries have been modelled as Lambertian on a 50% rough surface 2 .
Material VUV reflectivity Visible reflectivity metal 25% [45] 60% [45] reflector foils 0% [44] 93% [44]  We test our model in two different geometries, corresponding to two different experiments employing LArTPC detectors: subset-of-DUNE-like and SBND-like. The main parameters describing these two geometries are listed in Table III. In both geometries the PDs (PMT-like for the SBND-like geometry and X-ARAPUCA-like for the DUNE-like geometry) are distributed approximately evenly and with a PD located in the exact center of the PD-plane.  We simulate energy depositions at points within the active volume that cover the full phase space of distances and angles between the scintillation and the PDs. To study border effects and cover different regions of the detector, we divide the volume into concentric cylinders at different radial distances, d T , starting from the PD at the center of the PD-plane (Y − Z) outwards towards the corners of the field-cage, as illustrated in Fig. 3. Then at each d T , energy depositions are simulated at evenly spaced positions in the drift direction, X, covering the full drift length. In the SBND-like geometry, the simulated energy depositions are spaced in approximately 20 cm steps in the drift direction and 50 cm steps in d T . In the DUNE-like geometry, energy depositions are spaced in approximately 25 cm steps in the drift direction and 100 cm steps in d T . The step sizes in d T are driven by the distances between the PDs, while those in the drift direction have been roughly chosen to cover the entire detector volumes without requiring an excessive total number of points. Five different energy deposition locations are simulated for each X and d T pair, which define a cylindrical shell in three dimensions. The first 2 In Geant4: GLISUR with polish = 0.5 for the argon-metal boundaries and GroundFrontPainted with σ α = 0.8 rad for WLSfoil boundary. of these is chosen to be directly in front of a PD and the remaining four placed at increasing offsets of approximately 5 cm in d T to ensure a sufficiently large population of photons incident on the PDs from small angles. Following the above method, we have simulated approximately 2500 and 3500 energy deposition points in the SBND-like and DUNE-like geometries respectively.
Energy depositions generating 25 × 10 6 photons (∼1.04 GeV) are simulated in the SBND-like geometry and 100×10 6 photons (∼4.17 GeV) in the DUNE-like geometry. The larger number of photons in the DUNE-like case was chosen to improve the statistics in the number of photons incident on the PDs from the larger possible distances in this geometry.

IV. PREDICTING THE NUMBER OF DETECTED PHOTONS
In this section we develop a model capable of predicting the number of photons arriving at any given photon detector, based on the size of an energy deposit and the position where it occurred. We first focus on the VUV light that travels to the PDs directly from the point of emission, which we call the direct component. Starting from basic geometric considerations, we first calculate the solid angle subtended by a PD in an ideal infinite detector. We then account for the presence of Rayleigh scattering in the propagation by implementing a correction to the amount of light predicted geometrically. Finally, we add an extra correction to account for border effects present in real detectors of finite size.
The model developed to predict the number of VUV photons is then used as the first stage of a second model that predicts the number of photons arriving at PDs after being converted to visible wavelengths and reflected by a WLS-coated reflective detector cathode. We call this the reflected component of scintillation light. This model also employs corrections for transport effects and takes into account the finite size of the detector.
A. Direct VUV Light

Geometric Considerations
Scintillation light is emitted isotropically. This means that in an ideal case the number of photons arriving at a given PD could be calculated by simply estimating its geometric acceptance with respect to the scintillation point, i.e. the solid angle subtended by the sensitive window. This is a calculation that can be performed either analytically or with simple numerical integration depending on the shape of the PD. Such solutions exist, e.g. for disk [46] and rectangular [47] shapes, that cover most existing PD designs. The conclusions drawn in this work can be extrapolated to any other PD shape, once the calculation of the solid angle subtended from a point-like source by such a shape is known.
This approach works in an idealized case: in the absence of Rayleigh scattering and reflections by the detector materials (i.e. λ RS → ∞ and all materials 100% absorptive). In this scenario, the calculation becomes purely geometric and for any given energy deposition, ∆E, we can calculate the number of photons incident on a PD as, where the S γ (E ) is the scintillation yield of LAr for a given electric field, and Ω is the subtended solid angle. We also implement absorption effects due to contaminants: Q abs = e − d λ abs , where λ abs is the absorption length and d is the distance to the PD. The performance of Eq. 3 at predicting the number of photons can be seen in the top panel of Fig. 4. This shows a comparison between the number of photons hitting the PD windows predicted by Eq. 3 and the number obtained from a full Geant4 simulation, normalized to the sensitive-window area and the energy deposited. The pure-geometric calculation agrees with the full simulation within expected Poisson fluctuations. The gradient-colors of the circles represent the offset angle, θ defined in Section II E. It can be seen that more light is observed by the PDs from emission points closer and more on-axis, as expected. The red dashed line indicates a perfect 1/R 2 behavior. Even in this simple case, it becomes clear that it is necessary to account for the offset angle as that can change the prediction for a given distance by up to two orders of magnitude.

Corrections to the Geometric Approach
The basic solid angle approach breaks down when Rayleigh scattering is introduced into the simulation. The VUV scintillation photons in LAr undergo scattering during propagation with a characteristic length, λ RS , that is small compared to the size of current and future LArTPC experiments. This alters the path of the majority of the photons, and consequently the number of them arriving at the PDs. Once Rayleigh scattering is included in the Geant4 simulation the distribution of points in Fig. 4-top is altered to that shown in the bottom panel. The Rayleigh scattering significantly alters the amount of light observed in the PDs, and it is therefore essential to account for its effect during simulation. These effects strongly depend on both the distance, d, and the offset angle, θ, of the PD relative to the light emission point.
To build corrections for the effects of Rayleigh scattering, we calculate the ratio between the number of incident photons from Geant4 simulation, N hits , and the geometric estimation from Eq. 3 projected on cos(θ), N Ω /cos(θ). For simplicity, we split the phase space into 10 • -wide bins in θ. The discretization in θ introduces a systematic ef-  5. Relation between the number of Geant4 simulated hits on the PDs and the pure geometric estimation described by Eq. 3, in the SBND-like geometry. The error bars represent the standard deviation of the distribution within each angular bin. A strong dependency is clear in both distance and offset angle. At any angle, the dependency of the ratio with distance can be accurately described by a Gaisser-Hillas function, as illustrated by the dashed curves. To avoid large divergences at big offset angles we have found it more convenient to work with the projected solid angle. fect in our model: a more (less) sampled choice would result in a more (less) accurate correction. Our choice in this work is a trade off between accuracy and computational time. The resulting ratios, shown in Fig. 5, are smooth distributions as a function of distance and clearly separated between the different angular bins. This indicates that a parameterization in (d, θ) should include all the main dependencies and consequently be sufficient to predict the number of arriving photons. We find that the distributions shown in Fig. 5 can, for all angles, be accurately described using Gaisser-Hillas (GH) functions [48], as illustrated by the dashed curves: where N max is the maximum of the function located at a distance d max , and d 0 and Λ are parameters describing the width of the distribution. We implement the GH functions as the core of our numerical model to predict the scintillation light signals in large LArTPC detectors: (i) the number of incident photons on each PD is predicted by the solid angle that the aperture of the detector subtends, (ii) then the effect of Rayleigh scattering is accounted for via corrections to the geometric prediction.
Once we apply these corrections to Eq. 3, our model precisely predicts the number of incident photons on each PD as shown qualitatively in Fig. 4 (bottom panel, blue circles). A quantitative comparison is discussed in Sec-tionVI A 1.
In the limit d → 0, the effect of Rayleigh scattering should be negligible and the y-intercept in our corrections should correspond to the value cos(θ) (i.e. ∼1 for the on-axis case). The Gaisser-Hillas-like shape of the corrections suggests a behavior of the light such that for small distances from the PD, the probability of detecting scattered photons that would otherwise escape from the detectors is larger than the fraction that is lost due to the scattering. This situation continues for larger distances until a point at which it is reversed and more photons are lost than gained. Additionally, once taking into account the 1/cos(θ) factor in Fig. 5, we can see that PDs at large θ (that have a small geometric acceptance) present a higher relative probability to recover scattered photons compared with PDs located closer to on-axis. This behavior explains the significant tightening of the angular dependence in the points on Fig. 4 when Rayleigh scattering is included (bottom) compared to the ideal case (top), although the dependence remains strong. These effects also result in the detector size having an impact on the required correction curves: the greater the active volume in which photons can scatter, the greater the probability that these photons will end up feeding the signal. We account for this effect next.

Correcting for Detector Size: Border Effects
The dependency of our derived corrections on the detector size can be treated as a border effect. These borders (i.e. the cryostat walls and other detector components) not only delimit the active volume where photons can travel and scatter, but also consist of surfaces that can reflect or absorb them. These effects influence the amount of light observed in the PDs and, as a consequence, different sets of corrections may be required for different regions of the liquid argon volume.
To develop the corrections, we examine the behavior of the parameters of the Gaisser-Hillas functions as a function of the radial distance d T as defined in Section III. For simplicity, and taking advantage of the strong correlation between the d 0 and Λ parameters of the Gaisser-Hillas functions, we fix the value of d 0 absorbing all of the d T dependencies into the remaining three parameters. Figure 6 shows the results for the N max parameter (similar results are obtained for d max and Λ, and are shown in Appendix A). We observe a linear dependency in d T for all of the offset angle bins in both geometries under study. There is also a weak dependency in the slopes of the lines with θ, increasing for the more off-axis angles. We take these dependencies into account to accurately estimate the number of scintillation photons arriving at a PD for the entire active LAr volume. To this aim, we redefine the Gaisser-Hillas parameters in Eq. 4 as: where N max , d max and Λ are the values of the parameters in the center of the PD-plane (d T = 0 cm), and 1 , 2 and 3 are the slopes of the linear corrections for each parameter respectively. This new function is referred to as GH . To give an indication how these corrections affect the probability of photons arriving at PDs, Fig. 7 shows examples of correction curves for the two extreme cases, center (top) versus corner (bottom), for both the SBND-like and DUNE-like geometries.
Bringing all of the above effects together, the model we describe here is able to estimate the number of detected scintillation light photons using Eqs. 3 as: which depends only on the distance and angle between the emission point and the PD, and distance of the emission point from the center of the detector. Note that N γ = ∆E × S γ × Q abs × P (d, θ) × T (d, θ) using the notation from Eq. 1.

Basic Geometric Model
The number of photons arriving at the PDs in a LArTPC that has wavelength-shifting reflector foils on the cathode can also be predicted using a geometric approach. This requires expanding the model developed for the direct light VUV-only case described in Section IV A. The prediction of the wavelength-shifted and reflected visible light is inherently dependent on the specific detector geometry, because the distance between the reflective foils and the PDs becomes a key element of the model. Additionally, unlike for the VUV light, the wavelengthshifted light is much more likely to reflect from the borders of the detector and the field cage. Therefore, we construct the model for wavelength-shifted light using a realistic detector geometry from the start rather than using an idealized detector.
The visible light arrives at the PDs after being reemitted and possibly reflected by the WLS-coated reflector foils at the cathode of the detector. Therefore, we first calculate the number of VUV photons incident on the reflector foils using the solid angle that the entire cathode subtends, Ω c . This is corrected for the effects of Rayleigh scattering using Eq. 6. We then assume that these photons are re-emitted approximately isotropically after being wavelength-shifted and reflected, and that the region of the cathode in front of the scintillation in the drift direction will be the dominant source of the visi-ble photons. The central point of this region is referred to as the bright-spot, and is illustrated in Fig. 8 together with the other elements of the geometric model for the reflected light. The number of photons incident on each PD can then be calculated using the solid angle subtended by the PD aperture as viewed from the bright-spot, Ω P D . The geometric prediction for the number of visible photons arriving at the PDs can therefore be expressed as, where N γ,direct (Ω c , d c , θ c , d T ) is the prediction of the number of photons incident on the cathode using the direct VUV light model given by Eq. 6. The solid angle of the PD, Ω P D , is divided by 2π rather than 4π due to the presence of the highly reflective foils beneath the WLS.

Corrections for PD Position
The basic geometric model provides an initial approximation of the number of reflected photons incident on each PD. The assumption that the bright-spot region is dominant does not fully account for the distribution of the re-emitted wavelength-shifted photons across the whole surface of the reflective cathode. The approximation performs well for the PDs placed close to on-axis (at small θ c ) that see the majority of the light. However, it is a poorer approximation for the PDs located further off-axis where a larger fraction of the observed light originates from regions of the cathode opposite to the PD rather than the bright-spot. We therefore implement corrective factors to the basic model to account for these effects in an analogous way to the direct light model described in Section IV A. Because the corrections are developed in a realistic geometry, they also account for effects of reflections of the wavelength-shifted photons from the field-cage and the walls of the cryostat.
Similar to the method used for the direct light model, the required corrections are taken as the ratio between the number of incident photons on the PDs in Geant4 simulation and the prediction from the basic geometric model. For scintillation photons generated in the central region of a detector, the difference between the full Geant4 simulation and the model can be parameterized using only the distance between the scintillation and the bright-spot, d c , and the offset angle between the bright spot and the normal to the PD surface, θ c . Examples of the parameterized corrections are shown in Fig. 9 for the SBND-like and DUNE-like detector geometries. In both cases, the maximum offset angle as viewed from the bright-spot is defined by the size of the detector geometry. The error bars represent the standard deviation of the distribution within each angular bin. Its value can be affected significantly by the detector size, as can be seen comparing the DUNE-like case with the SBND-like one. This is caused by the large θ c angular bins describing larger regions of the detector volume where effects   arising from reflections on the detector walls can be significantly different. We note that the scintillation points with the highest variations, i.e. at high d c or θ c , account for a relatively small fraction of the total light observed.
To calculate the PD-position corrective factors we employ the mean of the N Geant4 /N Ω /cos(θ c ) distributions within each angular bin. Instead of using a fit to the cor-rective factors, we use linear interpolation in d c to find the exact correction for the prediction from the geometric model.

Correcting for Scintillation Position: Border Effects
In addition to the corrective factors accounting for the position of the PDs with respect to the point where the scintillation light was emitted, further corrections are required to account for the position where the light is created inside of the detector. Scintillation light created closer to the walls will be significantly affected by their proximity. The wavelength-shifted photons can be reflected off the walls, while the VUV light can be absorbed before it reaches the WLS-coated cathode plane.
We again account for these effects using parameterized corrective factors. We find that, similar to the PDposition based corrections, they depend on d c and θ c . An additional parameter is the position of the scintillation emission relative to the borders of the detector volume. Therefore, similar to the direct light model border corrections described in section IV A 3, sets of corrective factors are created at different distances, d T , from the center of the detector, as illustrated in Fig. 3. Then, during simulation, linear interpolation is used in both d c and d T for the required angular bin in θ c to calculate the exact corrective factor required.
Examples of sets of border effect corrections for the SBND-like detector geometry are shown in Fig. 10 for two cylinders defined by different values of d T . As before, the corrections are taken as the ratio between the amount of light seen in full simulation in Geant4 compared with the prediction from the geometric model. The equivalent corrections for the DUNE-like detector geometry are shown in Appendix B. The required corrective factors become significantly larger as d T increases and the scintillation is closer to the edges and corners of the detector volume. Additionally, the angular dependence becomes more significant and larger offset angles of the PDs, as viewed from the bright-spot, become geometrically possible. The border effects are much more significant for the light reflected by wavelength-shifting foils compared with the direct light and larger corrective factors are therefore required.
Bringing the above effects together, the number of incident reflected light photons on each PD can be expressed as, (8) where N Ω,ref lected is the geometric prediction given by Eq. 7 and A(d c , θ c , d T ) is the parameterized corrective factor accounting for PD position and border effects. This corrective factor depends only on the distance between the emission point and the bright-spot, the angle between the bright-spot and the PD, and the distance of the emission point from the center of the detector. As with the direct light, note that ∆E × S γ × Q abs × P (d c , θ c ) × T (d c , θ c ) using the notation from Eq. 1 for the reflected light.

V. PREDICTING THE PHOTON ARRIVAL TIME DISTRIBUTIONS
The model developed in the previous section addressed only the prediction of the number of photons arriving at each PD but not the distribution of their arrival times. As described in Sections II A and II B, the timing of the scintillation light is dominated by the double-exponential distribution caused by the de-excitation of the two argon molecular dimer states. In this section we describe a model to estimate the transport time t t (d, θ), see Eq. 2, that can affect the time distribution actually registered by the PDs. Analogous to Section IV we first develop a model for the direct light transport time and then use it as a starting point for a model describing light reflected off the cathode.

A. Direct Light Time Parameterization
The earliest arrival time of a photon on a particular PD can be predicted geometrically using the minimum distance that a photon must travel and the velocity of VUV light in LAr, shown in Fig. 1. A geometric calculation can provide the arrival time of the fastest possible photon, but does not account for other transport effects. A typical distribution of photon arrival times due to only transport effects can be seen in Fig. 11. The distribution shows a prompt component followed by a long diffuse tail.
We find that for essentially all combinations of emission point and PD the distributions are of a similar nature and can be approximated by a combination of Landau and Exponential functions using a simple fitting procedure. The resulting five parameters of the Landau + Exponential composite that describe a given time distribution are monotonic functions of the distance between the emission point and PD, provided we account for the incident angle. In the work described here we use two angular bins 3 : on-axis with θ ∈ [0 • , 45 • ], and off-axis with θ ∈ [45 • , 90 • ]. Figures 12 and 13 show the behavior of the model parameters for the two angle ranges in the two geometry cases: SBND-like (blue points) and DUNE-like (black points). The spread of the parameter values depends on the detector size: in a larger detector the signals are more scattered. For simplicity and because VUV photons are predominantly absorbed by all detector materials, we have neglected border effects in the model.
At larger distances the long diffuse tail of the arrival time distributions tends to disappear and the shape can be described using only a Landau distribution. We perform a quantitative comparison of the accuracy of the two approaches, as a function of the distance, using the relative difference of the χ 2 of both models. The result is shown in Fig. 14. In both of the geometry cases we find similar results: the Landau + Exponential model describes the shape of our signals more accurately, but at larger distances the two models perform similarly. The distance at which the two models become equivalent depends very slightly on the detector size, but for both geometries under study has a value around d = 400 cm. At longer distances the simpler Landau model can be used successfully instead of the Landau + Exponential one.
During simulation we construct and then sample the probability distribution of the VUV photon arrival times for each PD using the parameters of the Landau + Exponential composite function.

B. Reflected Light Time Parameterization
The transport time of photons arriving at the PDs as a result of a wavelength-shifting highly reflective layer on the cathode can be modelled using a similar approach. First, a geometric prediction of the transport time of earliest arriving photons is calculated. The fastest photons are most likely to travel along the path that mini- mizes the distance travelled at VUV wavelength, where the group velocity is slower. At visible wavelengths the photons propagate significantly faster due to the lower refractive index, as shown in Fig. 1. Figure 8 shows a diagram illustrating the most likely fastest path. The emitted VUV photons travel along the shortest path from the scintillation point to the cathode. There they are wavelength-shifted and re-emitted around the bright spot and take the shortest path to the PD. This simple model is able to predict the arrival time of the earliest photons.
The subsequent photons can be reflected from different regions of the wavelength-shifting foils and take very different paths to arrive at the PD. This results in a significantly broader distribution of their arrival times. We construct the model describing the visible photon arrival times at the PDs in three steps. We start by using the direct light Landau+Exponential model, described in Section V A, to estimate the arrival time distribution of the VUV photons at the bright-spot on the cathode. We then add the time needed for a visible photon to propagate between the bright-spot and the PD in a straight line. Finally, we apply a parameterized smearing to the result to account for the multitude of longer paths that can be taken. We use the following function for this smearing, where t s is the resulting smeared arrival time, t is the unsmeared arrival time, t f is the fastest possible arrival time calculated geometrically, τ is a smearing factor and x is a uniformly distributed random number between 0.5 and 1. This function keeps the earliest arrival times unchanged, but increasingly smears the photons arriving later. Additionally, a maximum time cut-off is applied to avoid an excessively long tail from the exponential distribution. We parameterize the smearing factor, τ , and the cutoff time, t max , in terms of the distance between the scintillation and the bright-spot, d c , and the offset angle, θ c , between the bright-spot and the PD, as shown in Fig. 8. The cut-off time is calculated as the time needed for 99.5% of Geant4 simulated photons to arrive. The τ parameter is determined by minimizing the difference between the smeared arrival time distribution calculated by the model and the distribution generated using Geant4. Unlike with the direct light transport time model, it is important to account for border effects such as reflections off the detector walls since they are highly reflective for visible photons. We use a similar approach to Section IV B, creating sets of smearing parameters at different distances from the center of the detector, d T . These sets can then be used to calculate the smearing parameters for any location in the detector using interpolation. An example set of the parameterized cut-off times and τ smearing factors is shown for the DUNE-like geometry in Fig.  15. An equivalent example for the SBND-like geometry is shown in Appendix B.
We observe that the cut-off times become larger with the size of the detector. This can intuitively be explained by the longer distances the photons need to travel before reaching the PDs, including many paths where they are reflected off the detector walls. The angular dependence of the cut-off time is relatively small, with a significant overlap between bins. The τ parameter is more dependent on the angle. This effect again grows with detector size and is much more prominent for the DUNE-like case.

VI. VALIDATION AND PERFORMANCE
To validate our model we test its performance against the results of a full simulation of the scintillation light in Geant4. For this test we use a sample of points generated in an analogous manner to that from Section III but shifted by several centimetres in random directions to test how the model works in the interpolated areas. We also compare the performance of our model with that of optical lookup libraries and give an example of applying the model to a realistic event.
A. Predicting the number of detected photons

Direct Light
The resolution obtained with the direct light semianalytic model as a function of d T is shown in Fig. 16 for the SBND-like and DUNE-like geometries. We obtain an unbiased estimation of the number of VUV photons arriving to our PDs in both geometries and for all values of d T . We also find the global resolution to be better than 10%, independent of d T . Figure 17 shows the performance as a function of the distance between the scintillation emission and the PD. The resolution worsens slightly with distance, ranging from 5-15% as we move from the closest to the farthest PDs. In each case, the performance is worst at the distances significant larger than the maximum drift distance (grey line) beyond which all PDs are off-axis. These PDs, however, are a minor contribution to the overall light signal of a physics event and do not significantly affect the overall resolution.

Reflected Light
The resolution obtained with the reflected light semianalytic model as a function of d T is shown in Fig. 18 for the SBND-like and DUNE-like geometries. The model performs well throughout the entire detector volume in both cases. It has a resolution better than 10% in the SBND-like geometry and better than 15% in the DUNElike geometry, with minimal bias in each case. For both geometries, the resolution is best in the central region of the detector, at small d T , where the effects of the borders are smallest. It then degrades slightly at larger d T as the border effects become more substantial and complex. The performance of the model in the DUNE-like certainty and spread in the corrective factors required, as seen in Fig. 9. Additional plots showing this effect can be found in Appendix B.
B. Predicting the photon arrival time distributions

Direct Light
The performance of the direct light model at predicting the time of the earliest arriving photon, t 0 , is shown in Table IV for the SBND-like and DUNE-like geometries. In both cases t 0 is predicted with a resolution better than 0.5 ns and with minimal bias. This resolution is smaller than the sampling of the PD electronics in current and upcoming LArTPC detectors, as described in Section II C. SBND-like detector geometry. The distribution of the photon arrival times is accurately predicted, except for a very slight tail-offset between the distribution from Geant4 and the model. This is due to the example θ lying at the extreme edge of the parameterized θ ∈ [0 • , 45 • ] angular bin. This offset could be reduced by increasing the number of angular bins used in the parameterization.

Reflected Light
The performance of the reflected light model at predicting time of the earliest arriving reflected photon, t 0 , is shown in Table IV for the SBND-like and DUNE-like geometries. In the SBND-like geometry, t 0 is predicted with a resolution better than 0.5 ns and without bias. In the DUNE-like geometry the performance is slightly worse, however the model still predicts t 0 with a resolution better than 1 ns and minimal bias. As with the direct light model, these numbers are well inside the timing resolution of the PD electronics in typical LArTPC detectors.
An example comparison between the reflected light photon transport time distributions predicted by the model and simulation in Geant4 is shown in Fig. 20 for the SBND-like detector geometry. The model accurately predicts the arrival time of the earliest photons and provides a reasonable approximation of their overall distribution. The model slightly underestimates in the first part of the tail of the distribution and overestimates towards the end of it. We found that this behavior most prominently affects off-axis PDs, which see substantially less light than those closer to the energy deposit, resulting in a relatively small overall impact.

C. Comparison with lookup libraries
An important consideration is how the performance of the model developed here compares to that of the lookup library method commonly used in neutrino LArTPCs. We perform this test for both the SBND-like and DUNElike detector geometries. To directly compare performance we generated dedicated lookup libraries with the same total number of photons used to train our model, see Table V for details. We used a uniform voxel size throughout the detectors and a uniform distribution of photons/voxel 4 . For completeness we also generated a "Hi-Res" version of the SBND-like lookup library to compare performance with a larger number of photons/voxel.
The results of the comparison for the direct, VUV light, model can be seen in Fig. 21 for the two geometries under study. We compare these plots to Fig. 17, where the performance for the semi-analytic model is shown (note the different axes). We find that our model behaves significantly better than the lookup libraries in terms of both bias and standard deviation, especially at larger distances. This is at least partially a result of under-sampling of the lookup libraries, as shown by the improved performance of the "Hi-Res" library in the SBND-like case. In the DUNE-like case the fluctuations are exacerbated by the fact that for distances larger than 450 cm the severe under-sampling in photons/voxel at the library generation stage causes the majority of predictions to be based on samples of less than 3 photons per voxel-PD pair. Additionally, at very short distances the lookup libraries suffer from a higher uncertainty due to the voxel size introducing discrete jumps in the predictions very close to the PDs. This second problem cannot be resolved by increasing the number of photons used per voxel, instead it requires reducing the size of the voxels or using a different approach altogether in this region. Figure 22 shows the performance of the generated lookup libraries in predicting the number of reflected light photons. Due to the nature of modelling the reflected light we use d T as the variable instead of distance from the PD. We compare these plots to Fig. 18, where the performance for the semi-analytic model is shown. We observe that in the SBND-like geometry the lookup library method performs comparably to the semi-analytic model, especially if the "Hi-Res" version is used, although a small under-prediction is observed in the regions of high d T . In the DUNE-like case the undersampling effects are so severe that the standard deviation of the lookup library prediction is much larger than for the semi-analytic method. The effect is again caused by many voxel-PD pairs where the prediction is made based on samples of a few photons. This could be mitigated by using a significantly higher number of photons/voxel to generate the lookup library.
Overall we find that the semi-analytic model performs significantly better than lookup libraries trained using the same number of photons.

D. Example Application to Realistic Events
We have shown that our model works well for predicting the number of photons and their arrival times from point-like energy depositions. In simulations of particle detectors we more often deal with "extended" objects such as tracks or showers. Our model can easily simulate these kinds of events using the paradigm used e.g in Geant4, where particle trajectories are composed of discrete energy depositions called steps. To simulate realistic particle events we can apply our model to each of these steps and combine the results to obtain the simulation of the full particle trajectory. An example of this approach is shown in Fig. 23, where we present the results of simulating the scintillation light originating from an anti-muon track decaying into a Michel positron inside the SBND-like geometry.
We also compare the prediction of the waveform observed by the PDs that our model makes to that of the full Geant4 simulation. We find excellent agreement for both the primary anti-muon scintillation peak and the secondary scintillation peak caused by the positron.

VII. CONCLUSIONS
We have presented a new method to predict the number of scintillation light photons incident on photon detectors and their arrival times that can be used for simulations in large liquid argon neutrino detectors. Two scenarios were considered: VUV scintillation light that propagates directly to the photon detectors and scintillation light that is reflected off a wavelength-shifter coated highly reflective cathode. In each case, the models start with a prediction from pure geometric considerations, then corrections are applied for photon transport and border effects. For the prediction of the direct VUV light, we obtain a resolution better than 10% in two different geometries: one SBND-like and one DUNE-like. For the reflected light, we obtain comparable performance in the smaller SBND-like detector and better than 15% resolution in the larger DUNE-like detector. In both scenarios, the prediction of the earliest photon arrival time provided by the models is within one nanosecond -better than the highest sampling used in liquid argon neutrino detectors to date. The method we propose is dramatically faster than the full Geant4 optical simulation and outperforms the currently used lookup library method when trained with the same number of fully simulated photons. It can be used in any large scale liquid argon detector with a simple tuning of the model parameters.  In this section we show additional tuning plots for the reflected light model. Figure 26 shows examples of the reflected light semi-analytic model border corrections at two different d T for the DUNE-like geometry. Figure 27 shows the reflected light transport time model cut-off times and τ parameters for the central region of the SBND-like geometry.