Space charge in liquid argon time-projection chambers: a review of analytical and numerical models, and mitigation method

The subject of space charge in ionisation detectors is reviewed, with particular attention to the case of liquid argon time projection chambers. Analytical and numerical description of the effects on the reconstructed coordinates along the drift and the transverse directions are presented. The cases of limited electron lifetime, of dual-phase detectors with ion feedback, and of detectors with small and comparable ratio between drift length and width are considered. Two design solutions that mitigates the effects are discussed.


Introduction
Ionisation detectors, based on gaseous or liquid media, are designed to collect the signal induced by the motion of electrons between polarisation electrodes. In comparison, the signal induced by ions is very low, because their drift velocity is lower by several orders of magnitude. However the low velocity also implies that the ions remain longer in the drift volume: in case of long drift distance or large irradiation, and in particular for liquid detectors, the charge density due to ions may be large enough to affect the electric field that drives the motion of the electrons. 1 Effects of space charge have been studied in calorimeters using liquid krypton [1] and liquid argon [2,3]. More recently, the relevance and the scope of this subject has widened with the development of large liquid argon imaging devices [7,8,9,10], where the drift volume may be long enough to become sensitive to the space charge due just to the exposure to cosmic rays at the surface of the Earth.
The purpose of this paper is to review the basic formalism describing space-charge effects in ionisation detectors. First, the simplest one-dimensional treatment will be reviewed, followed by a discussions of the effects of charge-yield dependence on the electric field, and of electron attachment, with the related presence of negative ions. Multi-dimensional cases are discussed in terms of boundary conditions on the side walls, and of aspect ratio (depth vs. width) of the drift volume. The case of dual-phase detector with injection of positive ions at the anode end of the drift volume is discussed. Two mitigation methods are considered, dealing respectively with the effects on the drift coordinate and on the transverse coordinates. Throughout the work, scaling laws are discussed, aiming at relations that could be used as guidelines for detectors operated in different conditions and with different geometry. Comparison between analytical approximations and numerical solutions are presented.

Review of basic assumptions and simplest case
This section will first summarise the results of [1], dealing with the simplest, one-dimensional equation for space charge and electric field. Next, it will consider aspect related to an ionisation detector used as a time projection chamber (TPC).

One-dimensional analysis
An ionising particle creates pairs of electrons and positive ions, but because the respective values of drift velocity differ by typically by 5 or 6 orders of magnitude, under steady conditions the density of the two charge carriers differ by the same amount. Therefore, in the fundamental approximation the density of charge is assumed to be due to positive ions alone, ρ + , which varies because of the injection of ions from ionising particles, and satisfies the continuity equation

∂ ρ ∂t
The stationary solution (∂ ρ + /∂t = 0) is of interest, and it is determined under the assumption of constant and uniform charge density injection rate K. The value of K depends on the flux of ionising particles crossing the detector and on the value of the electric field, which affects the initial recombination of electrons and ions. The space charge causes a non uniformity in the electric field, and therefore a space dependence of K, which is ignored here, and considered below in Sec. 2.3. The assumption of uniformity and stability of K is valid because the charge injection, due to e.g., cosmic rays, is effectively averaged over the time needed for the ions to drift from anode to cathode. The detector is taken as a parallel-plates ionisation chamber with gap L, operated with voltage V • and average electric field E • = V • /L directed along +x, with x = 0 (x = L) at the anode (cathode). The problem is reduced to one-dimension under the assumption that far from the side wall the border effects are negligible, so that the dependence on y and z may be ignored.
Under these assumptions, and using the mobility coefficient µ + , the continuity equation is solved as where v + x = µ + E x is used and the boundary condition ρ + (0) = 0 is applied, since the positive ions drift away from the anode and no accumulation of space charge is possible at x = 0. For x > 0, space charge is present reflecting the rate of charge density injection, accumulated over a time effectively equal to x/(µ + E x ).
The electric field satisfies the Gauss's law, which under these assumptions is written as and is solved directly as where E a is the value of the electric field at the anode, which is determined by the boundary E x dx = V • integrated from anode to cathode, while the adimensional parameter α is defined as [1] This parameter can be interpreted as α 2 being equal to the charge density injection rate K, multiplied by the gap length L and by the ions drift-time across the gap L/(µE • ), and divided by the surface charge density σ • = εE • at the electrodes, with the last two quantities computed for vanishing K. The ratio σ • /L ≡ ρ • is the natural unit for evaluating ρ + (x), which can be written as From Eq. 5, and as shown in Fig. 1, the field at the anode E a is always lower than E • , while the opposite holds for the field at the cathode E c . Correspondingly, the electric potential, defined as V (x) = − E x dx and integrated from the anode towards the cathode, deviates from the value −E • x, according to a smooth convex profile. The analytical expression for V (x) is given below in Eq. A.1, Appendix A.
As discussed in [1], a critical situation occurs for α ≥ 2, when the electric field vanishes at the anode, enhancing recombination between electrons and positive ions. 2 As an example, in a liquid argon detector operated at the surface of the Earth, the charge density injection rate due cosmic rays is approximately given by K = 2 × 10 −10 C m −3 s −1 . For L = 4 m, E • = 500 V/cm, 2 In the critical condition α = 2 the electric field varies linearly as E x = 2 E • (x/L) and the density of positive ions is uniform: For α > 2 the active region is reduced to a gap of length L = 2L/α, detached from the anode by The occurrence and the modality of a critical condition, predicted in [1], has been observed in calorimetric cells [3].
x/L under similar assumptions. For comparison, the curves corresponding to α = 1.6 are also shown.
It should be kept in mind that there is uncertainty in the value of the mobility of Ar + ions, with reported values in the range of 0.8 to 2.0×10 −7 m 2 s −1 V −1 [3,8,11,12,13], for values of temperature and pressure usual for particle detectors. Such uncertainties propagates directly to the value of the parameter α 2 , and to the application of the results presented in this paper.
Appendix A provides a set of relations valid in first-order expansion in the parameter α 2 , together with numerical approximations at higher order.

Time projection chambers
In a time-projection chamber, like liquid argon devices designed for neutrino detection [14], the effect of the space charge on the collection time of a charge deposited at a distance x from the anode is described by a time offset δt, given by where v e (x), v e • are the values of the drift velocity of electrons at the local electric field E x (x) and at the nominal value E • , respectively. Since the electron drift velocity depends monotonically on E x , the integral receives contribution of different sign from regions of small and large x.
to describe the non-linearity of the drift velocity of electrons, in the approximation of small values of α it is straightforward to find which naturally vanishes at x = 0 and also at x = L, although the latter property holds only to first order in α 2 . The maximum effect on the drift time occurs for x max L/ √ 3, and is equal to The effect on the reconstructed coordinate, which can be referred to as longitudinal distortion, is δ x = v e • δt, which in the same approximation has the largest value This equation shows that the longitudinal distortion is proportional to the charge density injection rate and to L 3 , E −2 • . In the example considered above (liquid argon TPC with α = 0.78, L = 4 m), and using γ 0.5 [15], the maximum longitudinal distortion is δ x max 7.7 cm.

Electric field variation and ionisation yield
The analytical description contained in Section 2 can be brought to a more realistic condition by introducing an x dependence in the change density injection, related to initial charge recombination, which changes the amount of free electrons and ions as a function of the electric field strength. To this purpose K is multiplied by R(E) and, with E = E x , Eq.s 2 and 5 are replaced by: Figure 3 shows examples of numerical solutions for E x (x)/E • . R(E x ) is taken from [16] and normalised to the value for E • = 500 V/cm: R(E x ) = 1.15/(1 + 72.9/E x ), with E x expressed in V/cm. Values α = 0.8 and 1.6 are considered, including recombination (continuous lines) and excluding it (dashed lines). The difference is negligible for the smaller value of α, and rather small for the larger. The effects of recombination would be more visible for larger values, and the threshold for critical density is increased to α 2.5.
Besides the effect on δ x, the x dependence of the field strength needs to be taken into account when the specific energy loss dE/dX along the trajectory of a charged track is extracted from the ionisation signal dQ/dX, for the purpose of particle identification. The local dependence of the charge yield affects the ionisation signal as dE ∝ dQ / R(E(x)). Furthermore, the segment on the track length dX = (dx 2 + dy 2 + dz 2 ) 0.5 requires a scale correction for the dx component, proportional to d(δ x)/dx. Using Eq. 9 and δ x = v e • δt, at first order in α 2 the correction varies between −0.08 α 2 cos 2 θ x at the the anode to +0.17 α 2 cos 2 θ x at the cafthode, where θ x is the angle between the track segment and the drift direction.
In the following sections, analytical results are obtained ignoring the dependence of the initial recombination on the electric field, while the effect is included in numerical computations.

Electron attachment and negative ions
Drifting electrons may be captured by electronegative impurities according to dρ e /dt = ρ e /τ e . The electron lifetime τ e is related to the attachment rate constant k e and to the density of electronegative impurities n i as τ e = 1/(k e n i ). The attachment rate constant depends on the type of impurity, and on the electron energy distribution, which is affected by the electric field, but this dependence has been observed to be rather small for field strength below 1000 V/cm [17].
The continuity equation for the average electron charge density, in steady conditions, can be written as where −K this time defines the injection of negative charge density, drifting towards x = 0 (with ρ e < 0, v e x < 0, K > 0). An average capture length λ e = |v e |τ e can be defined, and compared to the gap lenght L. The electron attachment is source of negative ions, which drift towards the anode with a mobility parameter similar to the one of negative ions [12], and contribute to the space charge present in the detector. The result of a numerical solution [18] to the continuity equations for ρ e , ρ + , ρ − and to the Gauss's low for E x is shown in Fig. 4 in the case of λ e (E • )/L = 2.58, α = 1.15, and assuming equal mobilities for negative and positive ions. The presence of the negative ions causes the minimum of the electric field to move away from the anode to a shallow minimum at x min /L 0.15 , where the field is 3% lower than at the anode; in comparison with the case of infinite lifetime, the field at the anode is increased from 0.77 E • to 0.83 E • , and the field at the cathode is barely changed from 1.38 E • to 1.37 E • . For an electron lifetime shorter by a factor 2 (i.e. λ e /L = 1.29), the corresponding value are x min /L 0.24, E min /E a = 0.94, An approximate analytical expression for E x (x), which matches the numerical evaluation at the level of 1%, is provided in appendix B. In the following sections, analytical results are obtained ignoring the dependence of the initial recombination on the electric field, while the effect is included in numerical computations.

Dual-phase detectors and feedback of positive ions
Dual-phase scintillation and ionisation detectors have been used for the study of rare processes, and large devices based on argon have been proposed for rare processes and neutrino experiments [19,20]. After drifting in the liquid, the electrons from primary ionisation are extracted into the vapour phase, where they can be accelerated in order to produce a light signal or to achieve charge amplification [21]. In the latter case, a significant fraction of positive ions from the multiplication process may be fed back to the liquid, contributing to the space charge in the drift volume.
Near the liquid-vapour interface, an extraction grid is designed to establish an electric field strength that facilitates transition of the electrons from the liquid to the vapour. At the interface, dielectric polarisation attracts charges in the vapour phase -of any sign -towards the interface surface, and repels them from the surface when they are in the liquid [22]. The extraction potential for electrons is designed to counteract this effect, together with any effective binding potential for conduction electrons in the liquid phase [23]. It has been argued that the discontinuity of the polarisation field at the interface may prevent positive ions from entering the liquid [19]. Opposite arguments have also been presented [24]. Lacking direct evidence, it is assumed here that positive ions can reach and cross the vapour to liquid interface. Polarisation and binding potential effects appear more likely to play a role in preventing negative ions, which may come from electron attachment, from leaving the liquid phase [25], as they do for electrons. That would be a minor effect under the assumption λ e L and, furthermore, a build-up of negative ions on the liquid-vapour interface would affect more directly the electron extraction from the liquid phase, rather than the electron drift in the liquid TPC.
The relative amount of ions feedback to the drift region is described by the product β = (g − 1) × f , where g is the gain on the electron signal and the factor f includes the collection of the positive ions from the amplification region, their transfer into the liquid phase, and the limited transparency of the extraction grid for positive charges. The value of β depends on the electric field at the grid, on the side  of the drift volume, which corresponds to the the field at the anode E a of single-phase detectors, and is affected by space charge.
In steady conditions, and neglecting electron capture for the moment, the flux of electrons leaving the liquid phase is equal to the total rate of primary ionisation in the detector, J e = KL, and the corresponding flux of positive ions crossing the extraction grid is The differential equation describing the space charge in steady conditions is the same as Eq. 2, but now the boundary condition includes J(0) and the solution is The differential equation describing the electric field can be be easily integrated also in this case, obtain- The effect of space charge are now controlled by the two adimensional parameters α and β . The reduction in E a as a function of α is significantly enhanced by the presence of feedback, as shown in Fig. 5. For β ≥ 1, the critical condition of vanishing E a is reached for values α < 1. Figure 6 shows the behaviour of E x /E • and ρ + (x)/ρ • for some values of the parameters α and β . While for β = 0 the space-charge density ρ + (x) increases with x, the trend may be inverted when ion feedback is present, as shown in these examples.

Mitigation technique n.1 : separation grid
At first order in α 2 , Eq. 4 can be written as , so that the difference in field strength between cathode and anode is approximately equal to This range of variation in field strength can be reduced by means of a third electrode, a grid placed at the coordinate x g that constrains V (x g ) to a suitable value V g . The grid generates a discontinuity in E x that can be exploited to increase the field strength at the anode and decrease it at the cathode. With a higher field on the anode side of the grid, drifting electrons will cross it, while a fraction of positive ions are captured, reducing further the effects of space charge in the region between the grid and the cathode. Figure 7 illustrates the value of E x (x)/E • for a given configuration of x g and V (x g ), and different conditions of space charge, corresponding to α equal to 0.8, 1.6 and 2. The grid reduces by at least a factor 2 the range between the highest and the lowest values of the field strength across the full gap. Figure 8 shows the corresponding distribution of ρ + (x)/(ρ • ). For x < x g , the electric field is described by Eq. 5 as in the case without grid, but now the boundary condition on the field at the anode E a is with δ g = −V g /V • and, as usual, V = −V • (0) at the cathode (anode). The boundary condition can be written as As discussed below, a convenient configuration is with x g /L = δ g , in which the grid restores at x g the voltage that would be obtained without charge injection. Then the boundary condition becomes showing that E a , together with E x (x) for 0 < x < x g , reproduces the solution for 0 < x < L obtained in the case without grid and with α replaced by α × x g /L. Therefore, the distortion to the values of the electric field are reduced by a factor equal to the square of x g , and the longitudinal distortion in a TPC device by the cube of x g .  As shown in appendix C, for x > x g , the general solution for the E x (x) is where E g− and E g+ are the field values at the grid, on the sides of x < x g , x > x g respectively, and are determined by the values of x g and δ g , together with those of E • and α. The flux of positive ions ρ + v + is reduced by the factor E g+ /E g− as it crosses the grid. The same factor applies for the change in ion drift velocity, so that derivative dE x /dx and the density ρ + are continuous across the grid, as shown in Fig.s 7 and 8.
The parameter x g and δ g can be chosen so that in the two regions defined by the grid, the ranges in the value of E x are approximately equal. This is approximately equivalent to aiming at the maximisation of the lowest value of E x across the full gap. As discussed in more detail in appendix C, choosing δ g = x g /L, i.e. restoring at x g the voltage that would be obtained without space charge and without grid, brings in proximity of the optimal configuration, and may result of practical implementation. The optimal value of x g /L is in the approximate range 0.6-0.7, with limited sensitivity to the exact value, apart from a preference for the lower (higher) values for α larger (smaller) than 1.
As shown in Fig. 7, the grid changes the range of E x − E • by a factor approximately equal to (x g /L) 2 0.5. The effect is larger when critical conditions are approached: for α = 2, the grid reduces the electric field variation from ±100% to about ±30%.
As a final remark, the separation grid is charged with negative surface charge density σ g = −ε(E g− −E g+ ) and subject to the electrostatic pressure p g = −ε(E g− 2 − E g+ 2 )/2. To first order in α 2 , the pressure is given by the expression p g −0.24 ε E 2 • α 2 , which is accurate within 5% for the examples shown in Fig.s 7 and 8. In practical situations it may be similar to the pressure on the anode (p a 0.5 ε E 2 • (1 − α 2 /3)), and is smaller than the pressure on the cathode (p c −0.5 ε E 2 • (1 + 2 α 2 /3)).

Side walls, field cage and transverse effects
Besides the effect on the component of the electric field driving the drift of electrons and ions discussed above, which can be referred to as longitudinal, in practical conditions the space charge can induce transverse distortion. The main reason for this effect is that liquid-argon time-projection chambers [7,8] have been operated with side walls equipped with field-cages designed for the operation without space charge, i.e. with a uniform gradient of the voltage V fc established by the field cage (dV fc /dx = −E • ). This field pattern does not match the field established far inside from the field cage, and a transverse component of the electric field arises close the side walls.
A constraint on the transverse components of the electric field E y , E z can be placed considering lineintegrals E ds computed along closed paths. For the component along y, consider the path of four straight segments shown in Fig. 9, which starts at a point (x, y, z) far from the field cage, (a) reaches the field-cage at (x, 0, z), (b) reaches the anode at (0, 0, z), (c) follows the anode to (0, y, z), and (d) closes the path to (x, y, z). The contribution from (c) is null, so that the opposite of the term in (a) is equal to the sum of the term in (b) and (d), and satisfies on the line-integral condition where the electric voltage at the field-cage is −E • x, the dependences of V (x, y, z) on y, z can be neglected because of the distance from the field cage, and V (x) is the electric voltage in the one-dimensional description of space-charge effects obtained from the integration of −E x (x) in Eq. 5. Since for 0 < x < L the absolute value of V (x) is smaller than E • x, the transverse electric field is negative (directed towards the field cage), and the drifting electron are focussed toward the center of the detector. The value of V (x) and and its approximation to first order in α 2 are provided in Appendix A. For α < 1.5, the absolute value of the integral is largest at x L/ √ 3 and is approximately equal to 0.064 α 2 E • L. For α approaching 2, the maximum is moved towards x L/2 with the value 0.25 E • L. The solution of a set of dependent differential equations describing the electric field and the motion of charge species can be numerically approximated using FEA software [18]. Computations have been made in three dimensions, or in two dimensions when the third coordinate is far from the field cage. Figure 10 shows an example of contours of equal values of E x and E y in a two dimensional computation.
If the transverse size of the detector is large compared to the anode to cathode distance L, the fundamental parameters that define the electric field are the same as in the one-dimensional case: E • , L, and α 2 , with the latter proportional to (L/E • ) 2 . Therefore, (E x − E • )/E • and E y /E • scale as α 2 , and L 2 /E 2 • , spatial distortions δ x, δ y in a TPC detector scale as α 2 L, and L 3 /E 2 • , when all quantities are computed on coordinates that scale as x/L and y/L. Numerical computations show that these scaling relations remain valid also for finite electron lifetime, if λ e ≡ |v e (E • )|τ e L. At the field cage, the longitudinal component E x (x, 0) is constrained to E • . For y > L/2 and x L/2, E x (x, y) is described by the one-dimensional approximation (Eq. 5) within about 1%. Near the electrodes, the one-dimensional approximation is valid for y > L to the same accuracy.
A finite electron lifetime, if significantly larger than the electron maximum drift time, does not alter the scenario, as shown in the comparison of Fig.s 10 and 11. The latter is obtained from a numerical calculation using λ e • /L = 2.58. As discussed above in Sec. 2.4, the presence of negative ions reduces the distortion in E x (the range −22% to +38% is reduced to −17% to +34% for α = 1.15, with the minimum of E x moved to x 0.15). The maximum values of |E y | is reduced by about 10%. The lateral extension of the region with significant transverse field is not significantly modified.
with the integration computed along the path followed by the drifting electrons. Under usual conditions, this transverse distortion collects same-sign contributions along the full drift path, and for an initial coordinate x > L/2, its value can be significantly larger than the longitudinal distortion of Eq. 11, because the latter is reduced by competing contributions of different sign from the regions x L/2 and x L/2, and also the reduced dependence of the electron drift velocity on the electric field (γ = (δ v e /v e )/(δ E x /E • ) 0.5 for E • 500 V/cm). The largest transverse distortion occurs for drift path originating at x = L and y = 0, where the numerical computation provides δ y max = δ y(L, 0, z) = 0.105 α 2 L, which is three times larger than the corresponding maximum longitudinal distortion (the coordinate z is assumed here to be far from the corresponding side walls). Drift paths are shown as black lines in Fig.s 10 and 11.
The transverse distortion δ y(x, y, z) reflects into a scale distortion along the y direction equal to d(δ y)/dy, which should be taken into account for measurements of specific signal yield dQ/dX, or when establishing intervals on the y coordinate. The largest distortion to the length scale dX occurs at the cathode and is approximately equal to −0.2 α 2 exp(−2 y/L) cos 2 θ y , where θ y is the angle between the direction of dX and the y axis.

Longitudinal electric field (V/cm)
Transverse electric field (V/cm) Fig. 11: Same as in Fig. 10, but now with electron attachment corresponding to τ e = 10 ms, and equal mobility for positive and negative ions.

Detector aspect ratio
If the transverse size of the drift volume (with widths W y , W z ) is not much larger than the drift gap L, the one-dimensional description of sections 2-4 needs to be revised. General features and numerical examples are presented in this section.
As discussed in section 5, with L W z and for the z coordinate far from the edges z = 0, W z , the transverse component of the electric field decreases as E y (x, y) ≈ E y (x, 0) × exp(−y/λ y (x)) as the distance y from the field cage increases, with λ y L/2. If W y is not much larger than L, there is not enough width to reach a negligible value of E y before approaching the centre of the drift volume y W y /2. Because of symmetry with the other half of the detector, E y still vanishes at y=0, but does it with a finite gradient The parameter κ vanishes at the electrodes (x = 0, x = L), and is a first order quantity in α 2 in the range 0 < x < L. Table 1 shows the result of a numerical computation of ∂ E y /∂ y for different detector aspect ratios, at coordinates x/L = 0.5, y/W y = 0.125, 0.25, 0.5, for α = 1.15, E • = 500 V/cm, τ e = 10 ms, and equal mobility of positive and negative ions. The dependence on z is ignored in the first three rows (W z L and L < z < W z − L), while in the last row W z = L and the gradient is shown at z = W z /2, where ∂ E z /∂ z = ∂ E y /∂ y.
Because of ∂ E y /∂ y, ∂ E z /∂ z, the Gauss' equation for the main component ∂ E x /∂ x used in Section 2 is modified. For a detector with W z = W y = W , near the axis y = z = W /2, both transverse component of the electric field contribute equally: Both ρ/ε and κ are positive and first order in α 2 , so that κ, which is due by the relative proximity of the linear field-cage, reduces the effects of space charge at the center of the detector. 4 The values shown in Table 1 should be compared to ρ/ε, which is equal to about 0.53 V cm −2 at the center of the drift volume.
Using the Gauss' equation, the different configurations shown in the table correspond to a reduction of ∂ E x /∂ x at the center of the drift volume, compared to the value obtained for L W x , W y , in the range of 10 to 65%.
The values of ∂ E y /∂ y shown in Table 1 can be scaled to different intensities taking into account that they are proportional to α 2 . Furthermore, for L = W y W z or L = W y = W z , and λ e L, the gap L is effectively the only parameter with the dimension of length, so that at the center of the detector: In the latter case, ∂ E y /∂ y is smaller because the amplitude of E y (x, y, z) is reduced by the constraints E y (x, y, 0) = E y (x, y,W z ) = 0 imposed by the field cage. However, the term ∂ E z /∂ z enters as well in Eq. 24. Naturally, the boundary conditions are more effective in reducing the effects of space charge when both lateral dimensions are comparable to the gap length. Analytical approximations to κ(x) are provided in Eq.s A.7, A.8 in Appendix A, together with approximations of the y dependence of E y .
The numerical computation of the field strength along the drift direction is illustrated in Fig. 12, for a drift volume with L = W y = W z = 6 m. The contours of equal values of E x are shown on the symmetry plane z = 3 m, and on the plane z = 0.5 m. The plot for z = 3 m can be compared to Fig. 11-left. The proximity with the field cage in both y and z reduces the range of E x on the detector symmetry axis to about 460-610 V/cm, a factor 0.6 smaller than in the case of large detector widths. At a z = 0.5 m the range of E x is reduced by an additional factor 0.5.
The narrow aspect ratio has a smaller impact on on the lateral distortion, since δ y is determined dominantly by the closest field cage. The comparison of Fig.s 12-left and 11 shows a 16% reduction in δ y(x, 0, z) at z = 3 m. Table 2 shows the comparison of different geometries, providing the extremes values for E x (on the detector axis, or far from the side walls), the maximum transverse component of the electric field |E y (x, 0, z)| max , and the maximum transverse shift of an electron path δ y (L, 0, z) max . The numerical values are provided showing explicitly the dependence on the scaling variables α, L, and can be applied to different configurations, with relative accuracy of some per cent for α < 1.5. The table also shows the 4 The continuity equation for the space-charge density is also affected by the parameter κ, but less directly than E x . Near the detector axis and considering positive ions, the transverse components of the current density satisfy ρ + v y + ρ + µ + κ (y − W y /2), and similarly for ρ + v z + . The gradient of ρ + vanishes along the y, z directions, so that the continuity equation becomes The second term on the right reduces the amount of space charge stored in the detector, but it contributes at order α 4 , with limited effects unless α approaches the critical value. The same conclusion holds if negative ions are taken into consideration.
x (m) x (m) result for a detector of equal sides (L = W y = W z ), including the case of positive ions feedback parameter β = 1, and the effect of the mitigation procedure discussed in Sec 7.
7 Mitigation technique n.2 : corrections to HV field-cage The transverse components E y , E z can be cancelled by means of setting the voltage gradient on the field-cage V fc (x) so that it reproduces the voltage profile V (x) corresponding to the one-dimensional description of the effects of space charge. In that condition, the field-cage is just as effective in removing border effects as the usual configuration of uniform gradient does for the case of α

1.
However the accuracy of the tuning V fc (x) is limited by different factors. Firstly, one has to rely on an a priori rather accurate knowledge of the ratio of charge density injection and ion mobility K/µ + , which includes the effect of ionisation yield and determines the value of α. A reasonable estimation of the electron lifetime τ e is also needed, if λ e is not much larger than L. Secondly, as discussed below in Sec. 8, convective motion related to thermal gradients and to fluid recirculation may affect significantly the distribution of space charge, introducing additional dependences on x, y and z, and preventing an accurate cancellation of the transverse components of the electric field.
Because of these inherent difficulties, an approximate correction to the voltage profile of the field cage is discussed here: besides the voltage imposed at the anode (V fc (0) = 0) and cathode (V fc (L) = −V • ), a third connection is provided to a resistive field cage at a coordinate x fc L/2, where the voltage is set a value V fc with absolute value lower than V • × (x fc /L) and closer to the voltage V (x fc ) observed at the same distance from the electrodes and far from the field cage. The voltage profile on the field cage remains linear, but with a change of slope at x = x fc . Differently from the usual configuration, where the drifting electrons are attracted inward along the entire drift gap, and in particular for x x fc , now the transverse component of the electric field is cancelled in the region where it was largest, and the remaining component at x x fc /2 and x (x fc + L)/2 have smaller amplitude.
The value of V fc can be considered as adjustable to the actual conditions of α, V • and, to some extent, to the effects of convective motion in which the TPC is operated, in particular if independent values of V fc can be chosen for four sides of the field-cage. Figure 13 shows the result of a numerical computation for a detector with L = 6 m, E • = 500 V/cm, and large width, like in Fig. 11. The (x, y) projection of the equipotential contours and of the drift paths are shown for the usual field cage, and for a correction applied at x = 3.6 m, where the voltage is set at −159 kV rather −175 kV. In this way, the transverse component of the electric field is locally cancelled. The maximum transverse distortion of a drift path across the full gap is reduced by a factor 2, which applies also at larger y values. The corresponding contours of equal E x and E y are shown in the bottom plots, which may be compared to those of Fig. 11.
For a narrow aspect ratio L W y or W z , as discussed in Section 6 the voltage profile established by the field cage affects the electric field across the entire detector volume. The correction to the field cage has the desirable effect of reducing the transverse components of the electric field, achieving a more uniform detector response in the plane y, z. On the other hand, it increases the range of the x-dependent effects, with a larger reduction in E x in the region of the anode and a larger increase at the cathode. The effect is shown in Fig. 14-top, which may be compared to Fig. 12. A correction at x = 3.5, with a reduction of the voltage from −175 kV to −161 kV, reduces the maximum transverse distortion by a factor 2, but the range in the longitudinal component at z = 3 m is increased by a factor 1.5.
A dual-phase detectors with feedback of positive ions requires larger corrections in order to cancel locally the transverse component of the electric field. Figure 14-bottom shows the case with β = 1, with the voltage at x = 3.5 raised to −141 kV. The comparison with Fig. 15 shows a reduction in the maximum transverse distortion by a factor 0.71, but an increase the range of variation of the longitudinal component of the electric field by a factor 1.4. The comparisons among different configurations are summarised in Table 2.
In the examples shown in Fig.s 13-14, the correction to the field cage is chosen aiming at a compensation of the transverse component of the electric field at x = x fc and y, z in proximity of the field cage. A overcompensation, with a larger correction to V (x fc ), may be used to invert locally the sign of the transverse component, obtaining a further reduction of the value δ y max for the full drift from anode to cathode. The limit to this option is posed the loss of drifting electrons from primary ionisation for x x fc and x (m) x (m)  sufficiently close to the field cage, which are driven outward and may be captured on the field cage without reaching the region x x fc /2, where they would be driven inward. If the active region at the anode starts at a distance ∆y from the field cage, a convenient design solution is to use a over-correction for which the maximum local outward displacement of drifting electrons matches the value of ∆y.
Considerations on the optimal values for position x fc and the voltage V (x fc ) are presented in Appendix D.
For detector with narrow aspect ratio, and in particular in case of feedback of positive ions, a more effective solution is the combined use of the mitigation technique discussed here together with the one presented in Sec. 4: a separation grid at x g 0.6 L, with V g = −(x g /L)V • , can be naturally coupled to a field-cage correction at the same coordinate and voltage. The range of E x and the value of δ y max would reduced factors approximately equal to 2 and 3, respectively, or larger, in case of major space-charge effects.
The fluid dynamics of large liquid-argon TPCs are studied with numerical evaluations [26]. The pattern of the liquid flow is affected by thermal gradients induced by heat transfer at the cryostat walls, by heat dissipated in electronics contained in the liquid -if present, and by liquid recirculation. The latter is performed for purification purposes, and contributes to fluid flow both directly and as an additional source of temperature non-uniformity. The value of the velocity field is typically predicted in the range of fractions of mm/s to several cm/s. For comparison, the drift velocity of positive ions in typical liquid argon devices is in the range of 5-10 mm/s. Therefore the space-charge density distribution and the pattern of longitudinal and transverse distortions may be altered in a significant way by fluid motion.
Large effects may be expected in proximity of elements that constrain the pattern of fluid motion, like side walls and possibly the electrode structures. Asymmetries and anomalous behaviour of transverse distortions δ y, δ z near the field cages have been reported [7,27], suggesting the relevance of convective motion, while others [8] have observed that at some distance from the field cage, the observed longitudinal effects δ x(x) are described in good approximation with the approach described in sections 2 and 6 without need to refer to fluid motion.
It is not clear yet how accurate can be a fluid dynamics model in predicting the distribution of space charge and the effects on drifting electrons, but studies are underway. Another desirable development would be a design of fluid recirculation, i.e. the geometry of inlets/outlets and the recirculation rate, that would take into account its influence on the distribution of ions, and minimise the related uncertainty in the prediction of space-charge effects.

Conclusions
The subject of space charge in large-size liquid argon TPC detectors has been reviewed, considering the effects on the longitudinal (drift time) and the transverse coordinates, and the implication on the measurement of the specific energy loss dE/dX. The subject is relevant for the important role that this detector technology is taking in present and future neutrino experiments.
Analytical description and numerical examples have been presented, underlying the dependence of the effects on detector size and operating conditions, and determining the adimensional parameters that drive the behaviour of the detector response. The potential enhancement of space-charge effects in dual-phase detectors with feedback of positive ions has been illustrated. Border effects, and the case of detectors with comparable longitudinal and transverse size have been discussed.
Two design solutions that mitigate the effects of space charge have been presented. In the simplest implementation, they can reduce by at least a factor 2 the longitudinal and transverse distortions, respectively.
The combination of the two solutions is a straightforward extension of the study presented here.

Appendices A Analytical expressions and approximations
In the one-dimensional, basic model in which the electric field is described by Eq. 5, the integral of E x (x) is given by: The value of E a is determined by V (L) = E • L. Figure 1 (and Fig. 5, for the case without ion feedback) shows the result of a numerical computation of E a /E • as a function of α. The following relations hold to first order in α 2 : anode, better than 0.01 (0.05) for α < 1.2 (1.6) , cathode, better than 0.01 (0.05) for α < 0.75 (1.12) , better than 0.0043 (0.033) for α < 1 (1.5) in full range of x.
All these approximations fail as α approaches 2. The maximum deviation between V (x) and −E • x is which occurs at x max between L/ √ 3 and L/2 for increasing values of α.
Approximations at higher order in α 2 have been obtained numerically:  At the field cage (y = 0) of a detector of large widths W y , W z L, and for z far from the field cage, the transverse component of the electric field is approximated by the expression with an accuracy of some per cent of the maximum value of |E y (x, 0)|.
The x dependence of the factor κ(x) describing ∂ E y /∂ y (x) at y = W y /2, z = W z /2 is given with an accuracy of a few per cent by the expressions: For W y L, W z L, L z (W z − L), x L/2 and y W y /2, the y dependence of E y (L/2, y) is approximated with an accuracy of a few per cent by E y (y) E y (0) × exp(−2.56 y/L − 0.32 y 2 /L 2 ) ≡ E y (0) × f L/2 (y) . (A.9) As discussed in Section 6, as y increases towards W y /2, the effect of the field cage at y = W y becomes more relevant, determining E y (x,W y /2) = 0. The combined effect of the two boundaries is well described by This equation is valid at the level of few per cent relative accuracy for L = 6 m, 0 ≤ x ≤ L and W y = 12 m, and better than 10% for W y = 6 m.

B Analytical approximation with finite electron lifetime in one dimension
In one dimension (i.e., far from the side walls) the set of equation describing the electric field and the charge density distributions for electron, positive and negative ions, including a finite electron lifetime, is: where ρ − , ρ e , v − and v e are negative, and ρ e is negligible when compared to ρ + and ρ − . The dependence of v e on E prevents a direct integration. An approximate solution can be found under the assumptions which brings to ρ e −K(L − x)/v e • , a reasonable approximation of the numerical solution shown in Fig. 4. The equations for ρ + , ρ − and E can then be integrated directly as where λ e • = v e • τ e . In the comparison with the numerical evaluation illustrated in Fig. 4, with α + = α − = 1.15 and λ e • /L = 2.58, the analytical approximation is found to be accurate to 1% in the values of E(x) and ρ + (x), and in the range (10-15)% in ρ − (x). For λ e • /L = 1.29, the accuracy is better than 2% and (10-30)% respectively. C Electric field with separation grid and optimal configuration at lowest order in α 2 As discussed in section 4, with a separation grid placed at x g with voltage V g = −δ g V • = −(x g /L)V • , the electric field for x ≤ x g is still given by equation 5: with the only difference that the value of the electric field at the anode E a is derived from the dependence of E a vs. α, shown in Fig. 1 and in Eq.s A.2, A.4, after replacing α with δ g α. Since at lowest order the effects of space charge are proportional to α 2 , in the region between the anode and the grid the distortion to the electric field is reduced by a factor δ g 2 . The value of δ g cannot be too small, though, in order to preserve favourable effects of the grid on the side x > x g .
In the region x ≥ x g , the boundary condition on the flux of positive ions is where K x g is the ions flux reaching the grid from the x < x g side, which cross the grid at a fraction equal to the ratio of electric field on the two sides. The continuity equation provides and the corresponding differential equation for E x (x) is directly integrated as given in Eq. 20.
The optimal grid position x g may be chosen so that in the corresponding case without grid, the field variation between x = 0 (anode ) and x = x g is equal to half of the total variation across the gap, so that the range of variation in E x in the two regions is approximately equal. At first order in α 2 , the position determined in this way is x g (1/ √ 2)(1 − α/16). Next, the value of δ g is chosen so that the range of E x is equal in the two region, namely: E a = E g+ , E g-= E c .
As an example, for α = 1.6, the parameters of the optimal configuration are shown in first line of Table C.1, together with the corresponding values of E x at the electrodes. The full range (E x|max − E x|min )/E • is 0.503, to be compared to 1.161 without the grid. The closeness between the values x g /L = 0.642 and δ g = 0.626 suggests to consider the configuration with equal values, which is shown in the third line of the table, and for which the range in E x /E • is increased by 0.052 only, and the lowest value of E x /E • is reduced by 0.034, which is also a small value, when compared to the difference of 0.270 obtained with the configuration without grid. Therefore the condition x g /L = δ g does not affect significantly the optimisation of the separation grid.
Additional examples that approach optimisation within the constrain x g /L = δ g are shown in other lines of Table C.1, together with a set of values obtained for x g /L = δ g = 0.7, which are used in Fig. 7. For α = 2, the ranges of E x /E • is 14% wider with x g /L equal to 0.7 rather than 0.6, but the the values E x|min /E • are similar, 0.63 and 0.65 respectively, so that in the end the effectiveness of the separation grid does not depend much on the choice of the value x g /L within the range 0.6 to 0.7. For a detector with wide aspect ratio L W y ,W z , and considering only positive ions, the voltage difference δV fc (x) between a point of drift coordinate x and transverse coordinates y, z far from the field cage and the corresponding voltage at the position x of a linear field-cage is given by the voltage distortion δV (x) in the one-dimensional description of the effects of space charge. Near the field cage, the voltage difference generates transverse components of the electric field according to Eq. 21, and the most effective position x fc for a correction to the voltage profile of the field cage is at the maximum of δV fc (x).
For α < 1.5, δV fc (x) = V (x) + E • x is approximated using Eq.s A.2, and the preferred coordinate is x fc L/ √ 3. A correction voltage V fc = −(E • L/ √ 3)(1 − α 2 /9) is applied to the field cage at that point in order to cancel δV fc locally and reduce it in the full range 0 < x < L. Compared to a field cage without correction, the current flowing from the anode to x fc is reduced by a fraction equal to α 2 /9, the current flowing from x fc to the cathode is increased by a fraction equal to α 2 /[9( √ 3 − 1)], with the difference of about 0.26 α 2 being provided at the intermediate connection.
Once the correction is applied, δV fc (x) takes a different shape and two local maxima occur at x = (1/3)L and x 0.80 L. The values of the maxima are equal to α 2 E • L/81 and about α 2 E • L/57 respectively, much smaller than α 2 E • L/16 obtained without correction to the field cage.
For larger values of α, without applying the correction, the maximum of δV fc (x) moves towards x = L/2. For instance, for α = 1.5, the preferred value is x fc 0.54 L, with the maximum of δV fc (x) still approximated by α 2 E • L/16 within about 1%.
When negative ions from electron capture are present, smaller values of δV fc (x) are found for the same value of α. Computation performed with τ e = 10 ms, L = 6 m, E • = 500 V/cm and α = 1.15, have shown that, without the correction, δV fc (x) is about 10% smaller than in in the case with infinite electron lifetime, with the coordinate at the maximum of δV fc still well approximated by x L/ √ 3.