Mathematical Model of Four-Line Probe to Determine Conductive Properties of Thin-Film Battery Electrodes

Measurement of the electronic conductivity of porous thin-film battery electrodes poses significant challenges, particularly when the film is attached to a metallic current collector. We have developed a micro-four-line probe and testing procedure that overcomes many of these difficulties while relying on principles similar to commonly used four-point probes. This work describes a mathematical model that enables rapid inversion of the data collected by such experiments to compute two properties: bulk electronic conductivity of the film and contact resistance with the current collector. The model accounts for variable probe and sample geometry and variable resistance between the probe and the sample. Results from 2D and 3D models are presented. The full 3D model combines a Fourier series with the boundary element method to generate a solution that requires significantly less computational cost than a corresponding finite element solution for the same level of accuracy. The model confirms that the ideal probe line spacing is close to the value of the electrode film thickness. © The Author(s) 2015. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0571510jes] All rights reserved.

Transportation, communication, and mobile electronics are just a few of the many industries that are relying increasingly on battery technology. Lithium-ion batteries that employ particle-based thin-film electrodes are a critical part of this technological advance. However, such electrodes can exhibit conductivity limitations, including significant spatial variations due to particle and composition inhomogeneities, which can cause so-called hot and cold spots, undesired sidereactions, and an overall decrease in cell performance and safety. [1][2][3] One key to resolving such manufacturing and materials problems is accurate measurement of relevant properties. The general approach for conductive properties is to perturb the sample with an imposed current and then measure the electric potential on two or more points on the surface. However, it is surprisingly difficult to accurately measure the electrical conductivity and contact resistance (between electrode and current collector) of these films. [4][5][6][7][8][9] This is because the contact resistance between probe and sample is high relative to the resistance of the sample, the sample is mechanically fragile, and the conductivity of the electrode film is confounded by the presence of a highly conductive adjacent metal layer (the current collector). The electrode film can be delaminated from the current collector in order to eliminate the last problem, 4 though it is more desirable to measure electrodes non-destructively as is proposed in this work.
Conductive property measurements.-Multiple methods have been used in the past to measure the conductivity or resistivity of thin films or layered materials. The four-point-probe method has been in use for many decades, and is particularly useful in measuring singlelayer conductive materials. On the other hand, with point probes it is difficult to control the pressure imparted to the electrode material during the measurement, which can change conductivity significantly, 4,10 Furthermore, with particle-based samples it can be difficult to apply the appropriate amount of point force, as too little force leads to poor probe contact and too much force can penetrate and damage the sample.
The van der Pauw method is another common way to determine the resistivity of a sheet regardless of its shape, as long as a few critical parameters are met: uniform resistivity, uniform thickness, no isolated surface holes, and the sensors being located on the circumference of the shape. 11 Electrical impedance tomography is a method of using multiple current-imposing probes and multiple voltage sensors across a volume. It is proposed for use in medicine as a non-invasive technique to determine tissue characteristics. Coupling an associated numerical z E-mail: dean_wheeler@byu.edu model with the experimental results enables one to map out the conductivity of the tissue sample. Accurate results depend on creating an accurate model, which is a difficult problem that requires a good deal of knowledge about the domain being sampled. 12 For example, if the sensor contact resistances are not accounted for in the model, or if the geometry and positioning of the electrodes are not accurate, the model will provide a flawed inversion and subsequently predict incorrect distribution of tissue properties.
Electrical resistance tomography is similar to impedance tomography, but with a focus on prospecting, mining, and other geologic applications. 13 Electrodes are placed in boreholes in order to detect underground objects by the resulting potential perturbations, again requiring numerical inversion techniques. A related imaging method has also been reported as being useful for process control in industrial equipment such as hydrocyclones. 14 Four-line probe.-The overarching idea of the above techniques is that imposed currents create variable surface potentials, which can be measured and used with a model to regress the resistance/conductance properties in the interior of the sample volume. This paper presents the mathematical development of a model that can be used to interpret experiments based on a four-line probe. The new probe extends the idea of a four-point probe and overcomes the aforementioned challenges of pressure control and measurement of a thin film still attached to a current collector. Particular attention is paid to developing a method that is robust and sufficiently sensitive, given typical noise levels, so that local values of conductivity can be determined for thin-film battery electrodes. In an accompanying paper 2 we show experimental results for a micro-four-line probe, constructed using semiconductor fabrication techniques, and test fixture that have been used to probe spatial variations in electronic conductivity and contact resistance for commercial-grade lithium-ion battery electrodes. The model presented in this paper is the first step in that effort. Fig. 1 shows the essential difference between a four-point measurement and a four-line measurement. Both accomplish the basic function of a four-sensor probe in mitigating measurement effects from probe contact resistance. However, the four-line probe also provides greater contact area between probe and sample while still maintaining a small spacing between the contacts (i.e. lines). Furthermore, the lines for our device are embedded in a flat insulating substrate, enabling pressure on the sample to be controlled and measured during the electrical measurement. 2 Conduction models.-A mathematical model is needed to extract a film conductivity from a four-sensor measurement. As shown in Fig. 1, a current I and a potential difference V are measured. Ohm's law leads to an apparent resistance R = V /I. This resistance in turn In this depiction, current I is injected and removed using the two outer contacts, while potential difference V is measured using the two inner contacts.
is proportional to the resistivity ρ of the sample (i.e. inversely proportional to the conductivity σ), with the constant of proportionality being a shape factor S. To determine the shape factor, Laplace's equation ∇ 2 φ = 0 for local potential φ must be solved for the particular geometry and boundary conditions used.
Traditional four-point-probe analysis requires some simplifying assumptions: infinitesimal area of contact, symmetric probe spacing, a single-layer film of infinite planar extent, and isotropic electronic conductivity. 15 The assumption of infinitesimal area of contact allows the current injection and removal (the two outer contacts) to be treated as point sources and the voltage sensing (the two inner contacts) to likewise be treated as discrete points. The absence of a current collector means the lower boundary can be treated as insulating. The remaining assumptions likewise facilitate a simpler mathematical solution. Under these assumptions and for a sample film that can be approximated as infinitesimally thick, the mathematical solution of Laplace's equation can be performed conveniently using the bipolar cylindrical coordinate system, 16 leading to shape factor S 4PP = ln 2 (πT ) −1 , where T is film thickness. Conversely, if the sample is approximated as a semi-infinite medium (extremely thick), a solution using the bispherical coordinate system 16 leads to S 4PP = (2πP) −1 , where P is the center-to-center distance between adjacent contact points. Films whose thickness is similar in magnitude to the probe spacing (T ≈ P) require a more complex solution.
For the four-line-probe model we assume an isotropic and uniform conductivity σ for the region of interest, i.e. the vicinity of the probe, even though conductivity may vary on larger length scales. It is likely that typical battery films exhibit some anisotropy in conductivity due to nonspherical particles and directional fabrication steps, namely coating, drying, and calendering. 17 Treatment of this is beyond the scope of the present work, but will be addressed in future work.
At this stage we seek a robust solution for multiple possible probe and sample geometries and for boundary conditions realistic for thinfilm battery electrodes. Typical Li-ion battery electrodes have film thickness on the order of 30 − 100 μm. The width of our probe is on a similar length scale, meaning the geometry is not amenable to simple analytical solutions. While one can use the finite element method (FEM) to generate a solution for such cases, it is generally quite computationally expensive to simulate the full 3D geometry. We wish for a more rapid solution so that the inversion procedure can be completed in real time using a desktop computer during sequential sample measurements.
Additionally, it is of interest to measure the area-specific contact resistance R c between the current collector and the electrode film. This adds an additional unknown and thus two independent experiments are required in order to reliably determine both σ and R c for the sample, as shown in Fig. 2. The first experiment (Fig. 2a) is akin to a traditional four-point-probe measurement ( Fig. 1) with the current largely flowing tangentially to the sample surface. With a conventional sample, this would simply measure σ, but the measurement is skewed in the case of intact battery electrode films by a shunt current that flows through the attached current collector. In fact, for traditional four-point probes in which P T , the majority of the current between the outer contacts would pass through the current collector, bypassing the region near the inner contacts and making conventional shape factors inaccurate for battery work. (For further details on this topic, see Fig. 5). The additional orthogonal experiment in our design ( Fig. 2b) causes the current to pass from the outer lines, which are electrically joined in this case, to the current collector, which is grounded. This second experiment provides a more direct measurement of the current collector contact resistance, though in general both experiments depend on both σ and R c and a simultaneous regression is required.
The remainder of the paper explains the development of two versions of the four-line-probe model: a two-dimensional (2D) model with simpler boundary conditions (BCs) and a three-dimensional (3D) model with more complicated BCs. Even though the 3D model is more realistic, it is more difficult to program and so the 2D model is valuable for testing and validation. Finally, the two models are compared in terms of resulting shape factors and potential distributions. The reader is referred to the accompanying experimental paper 1 for a description of a fabricated micro-four-line probe and measurements of the mmscale spatial distribution of conductivity and contact resistance for commercial-grade Li-ion electrodes, generated with the help of the 3D model.

Two-Dimensional Model
Boundary conditions.-The mathematical formulation of the 2D four-line-probe model begins with a solution to potential φ for the geometry and BCs given in Fig. 3. The length of the contact lines, L, is generally quite large relative to the dimensions shown in Fig. 3, justifying the use of a 2D cross-sectional model. A line of symmetry at x = 0 is employed to decrease from four to two the necessary number of contact lines included in the model. An insulating boundary occurs at x = W x , which is chosen to be large enough to approximate semi-infinite behavior. The Robin-type boundary condition at y = 0, R 0 T dφ/dy − φ = 0, reflects how some current penetrates into the current collector during the measurement. The metallic current collector is assumed in the model to be at a constant potential (ground) as any potential variations inside it can be shown to be negligible because it is so much more conductive than a battery electrode film. R 0 is a dimensionless surface impedance at the current collector, and reflects the relative resistances of the surface and the bulk film: (In this paper, sans serif font is used for dimensionless variables such as R 0 .) A mixed-Neumann BC is used on the top edge (y = T ). The inner line spans p 1 ≤ x ≤ p 2 and the outer line spans p 3 ≤ x ≤ p 4 . Away from these contact lines the surface is insulating. We assume a relatively large interfacial resistance between the contact lines and the sample (this assumption is relaxed for the 3D model).
For the outer line this assumption produces a uniform gradient or current density normal to the boundary. In effect, the large interfacial resistance equalizes the delivery of current through the contact surface. The value of boundary gradient G is immaterial to the resulting shape factors-the currents and potentials are each proportional to G. For the inner line the large surface resistance results in an insulating boundary. The floating voltage of the inner line, V inner , will be the average of φ across the width of the contacting the line.
Series solution.-The analytical solution to Laplace's equation for the 2D model is produced by standard techniques, 18 and so we only summarize the procedure here. The steps are (1) assume a solution by separation of variables; (2) satisfy the homogeneous BCs, leading to particular eigenvalues and eigenfunctions; (3) satisfy the final mixed/inhomogeneous BC at y = T by a Fourier series and the orthogonality of the eigenfunctions. Several components of the 2D solution are designed to be reused in the subsequent 3D model. The solution can be represented as follows.
where n = 0, 1, . . . , N and N is a sufficiently large integer to ensure convergence of the solution. The first component in the series solution is eigenfunction Y : This function is scaled so as to generate unit gradient on the top edge, (dY/dy) y=T = 1. Next, the eigenvalues λ n take different values depending on whether we are solving for the tangential (superscript t) or orthogonal (superscript o) BCs: There is a similar distinction for eigenfunctions X n : where δ n,0 is the Kronecker delta and serves to make X o 0 properly normalized so that the n = 0 term can be directly included with the others in the sum. Fourier coefficients C n for the tangential and orthogonal experiments are Arbitrary parameter G (Fig. 3) is built into these C n values. Functions Shape factors.-In order to determine conductivity σ and contact resistance R c , the model is compared to experimental results by means of shape factors for the tangential and orthogonal experiments. Each shape factor S is defined in terms of an apparent resistivity ρ that is a ratio of V inner (relative to ground) and I : The factors of 2 used above are designed for experimental convenience, and make no difference in the method as long as they are consistently applied in the shape factors and experimental measurements. As formulated here each S is dimensionless and a function of the geometry (T , p 1 , p 2 , p 3 , p 4 ) and ratio R 0 , but not of σ explicitly. The basic procedure is to measure the ρ t and ρ o quantities experimentally and then to solve for σ and R 0 by iteration: for each trial R 0 value one uses the model to compute S t and S o and iteration continues until Eqs. 11 and 12 are mutually satisfied. The final value of R c can then be determined from Eq. 1. The fact that shape factors must be computed repeatedly in the inversion procedure is part of our motivation to make their calculation as efficient as possible.
To compute shape factors from the model we must generate the ratios on the right side of Eqs. 11 and 12. Total current I for either the tangential or orthogonal case is obtained by summing the normal current flux across the outer-line contact region ( p 3 ≤ x ≤ p 4 , y = T ): V inner is obtained by averaging the solution to φ across the inner line region ( p 1 ≤ x ≤ p 2 , y = T ): Substituting Eqs. 13 and 14 into Eqs. 11 and 12 yields model expressions for S t and S o that contain R 0 and geometric parameters, but not σ or G.
In the limit that T → ∞ the 2D model for the tangential experiment can be solved analytically using the bipolar cylindrical coordinate system. The resulting shape factor is independent of the R 0 value: p + p p − p dp dp [15] This has a straightforward analytical solution, though it is algebraically lengthy and so is not presented here.
In a similar vein, it is possible to take the discrete sum of Eq. 2 and turn it into an integral by taking the limit W x → ∞. 18 This corresponds to a true semi-infinite boundary in x, rather than using a large W x value that imitates the same. However, this has little practical effect provided W x and N were sufficiently large in the first place (see Eqs. 44 and 46 below). Quadrature-type numerical evaluation of the resulting complicated integral may afford some computational efficiency over the discrete Fourier sum, but is not done in this work.

Three-Dimensional Model
Boundary conditions.-The 2D model solution can be computed extremely rapidly. However, its two drawbacks are neglect of the finite-length of the contact lines and the assumption of a large interfacial resistance between the probe and the sample. If the contact lines were infinitesimally wide ( p 1 = p 2 and p 3 = p 4 ) then the amount of interfacial resistance on the top surface would not affect the shape factor. However, for finite-width contacts, this resistance can make a small difference, even for a four-contact design that is meant to correct for this. Therefore the 3D model is an attempt to remedy these moderate drawbacks, at the cost of increased complexity and computational cost.
The geometry and BCs for the 3D model are given in Fig. 4. The end projection (bottom of Fig. 4) is similar to the 2D scheme in Fig. 3 and needs no further comment. The addition of the third dimension means that the domain extends to z = W z , some distance past the end of the lines at z = L. Likewise, the x y plane at z = 0 is a plane of symmetry, meaning that L is properly understood as half of the total (experimental) length of the contacts. In fact, the use of two planes of symmetry (x = 0 and z = 0) means the 3D model only needs to include one quadrant of the full probe/sample domain. The top surface of the 3D model contains a mixed BC. As with the 2D model, the surface is insulating where there is not a contact line: On the outer contact line ( p 3 ≤ x ≤ p 4 , y = T , 0 ≤ z ≤ L) a Robin-type BC is used: [17] where V outer is a fixed non-zero quantity. Because this is the only inhomogeneous boundary, all potentials and currents in the system scale proportionally to V outer , making it an arbitrary quantity that can be given a convenient value when calculating shape factors. R T is the dimensionless contact impedance between the probe and sample (defined in the same manner as R 0 in Eq. 1). It controls the distribution of local current flux across the contact area: if R T is large the injected current will be more uniform in the manner of the 2D model; if R T is small the current will concentrate more at the edges of the contact area and φ across the boundary will be more uniform, in the manner of a primary current distribution.
The BC for the inner contact line ( [18] where the second equality indicates that V inner is a floating potential that is determined on-the-fly by averaging φ over the inner contact area. This relationship is equivalent to stating that the total current for the inner contact is zero, though local currents can be nonzero. Furthermore, this means Eq. 18 is not a conventional Dirichlet, Neumann, or Robin boundary, although it is linear and homogeneous. The use of spatially constant V outer and V inner in the model is based on the reasonable assumption that the probe lines are internally highly conducting so that each operates at a single potential. Similarly, the metallic current collector is assumed to be at a constant zero potential for purposes of the model.
Series solution.-Initially the solution procedure to the 3D case follows the same steps as the 2D case, namely separation of variables and satisfying the homogeneous BCs first. 18 This leads to a doubleseries solution: φ(x, y, z) = n m C nm X n (x) Z m (z) Y (γ nm , y) [19] where n = 0, 1, . . . , N and m = 0, 1, . . . , M. N and M are sufficiently large integers to ensure satisfactory convergence of the solution. X n and λ n remain the same from the 2D solution. Function Y is still defined by Eq. 3 but uses a combined eigenvalue: where the eigenvalues for the z-direction follow the pattern for previously defined λ o n : μ m = πW −1 z m [21] Likewise, eigenfunctions Z m follow the pattern for X o n : FunctionZ m (used below) also follows the pattern ofX o n (Eq. 10). In contrast to C nm , X n , and γ nm in Eq. 19, values for μ m and Z m are independent of whether a tangential or orthogonal case is being computed.
Eq. 19 is guaranteed to satisfy all of the homogeneous boundaries (every surface except y = T ) for an arbitrary choice of Fourier coefficients C nm . To satisfy the inhomogeneous remaining BC at y = T , the proper choice of C nm values must be made. If the BC at y = T exactly specified the normal potential gradient, φ , at each point on the surface, then one could use an integral transform as was used for the 2D model, namely Unfortunately this type of analytical procedure for determining the unknown coefficients does not work for the mixed-type BC used in the 3D model. However, rather than abandon the series solution entirely, we attempt a numerical calculation of C nm by the boundary element method (BEM). 19 Our hybrid analytical/numerical method differs from, but has a similar spirit to, other semi-analytical methods that have been used to solve challenging boundary value problems. 20,21 Hybrid boundary element method.-With the BEM, the outer surface (in this case only y = T ) is divided into elements or nodes indicated by index i or i . Each element can be a source or sink of current, given by strength c i , which in general exerts an influence on the potential of all boundary elements in the system. The spatiallyaveraged BC equation for element i is put in the form where A ii is the so-called influence coefficient between elements i and i . The objective is to adjust the c i values until the BC equation at each element i is satisfied in an average sense. Because it is formulated as a matrix inversion problem, Eq. 24 is straightforward to solve. Therefore the principal task with the hybrid model is to relate the parameters in Eq. 24 to those in Eq. 19.
In this implementation the boundary elements are rectangular regions located in the vicinity of the inner and outer contact lines. Away from this region the surface is insulating and c i values would be essentially zero, so it is not necessary to place boundary elements over the whole y = T surface. Elements are not equally sized. Smaller elements are placed where potential variations are large, near the edges of the inner and outer contact areas. It was found sufficient for the calculations performed in this work to span the x direction with 14 divisions and the z direction with 9 divisions, making 14×9 = 126 total elements. Of those elements, 40 corresponded to the inner contact area, 40 to the outer contact area, and 46 to the nearby surrounding insulating surface. The region bounded by each element i is designated by In order to populate the element-averaged BC equations embedded in Eq. 24, the element-averaged potentialφ i , its normal gradientφ i = ∂φ i /∂ y, or both must be determined. These are expressed in terms of averages over the eigenfunctions in the x and z directions (compare to Eq. 19):φ where in Eq. 26 we take advantage of the fact that the derivative of Y is exactly 1 at y = T . H in and J im are constants that average sine or cosine eigenfunctions across a boundary element: In addition we define the following terms that when substituted into Eq. 25 allow us to evaluate V inner (see Eq. 18 and Eq. 36): Note that H terms differ according to the tangential and orthogonal BC cases, whereas J terms do not.
A key step of this hybrid method is to expand Fourier coefficient C nm in terms of boundary element strengths c i and averaged eigenfunctions, in a manner suggested by Eq. 23: If the Fourier series were not truncated (i.e. if N , M → ∞), then the solution to c i would represent exactly the total current flowing through element i. For a finite number of Fourier terms c i is still defined by Eq. 31, and may be considered a close approximation of the total element current.
One can now assemble the terms needed in the BEM Eq. 24. Eq. 31 substituted into Eqs. 25 and 26 gives quantities needed for each element's BC, which will be of insulating, outer-line, or innerline type. By comparison with Eq. 24, A ii and B i are then determined. For element i on an insulating region (Eq. 16) the result is For element i on an outer line (Eq. 17) the result is For element i on an inner line (Eq. 18) the result is While the equation for each A ii requires an expensive double sum, a judicious arrangement of the order of the calculations and reuse of terms can make the procedure efficient and manageable. Element strengths c i are then determined by a matrix solution procedure, such as Gaussian elimination. Because the size of influence matrix A in the hybrid method is much smaller than it would be for a standard BEM or FEM solution, this step is practically trivial.
Inner and outer shape factors.-Shape factors for the 3D model follow the pattern of the 2D model (Eqs. 11 and 12), except that they are also functions of R T . For the inner line they are Because there are now three unknowns (R T , R 0 , and σ), at least one additional piece of information must be collected from the same tangential and orthogonal experiments. We therefore devised two additional shape factors S * based on the outer-line voltage: The additional shape factors (denoted with superscript *) are particularly sensitive to R T and are therefore complementary to the original inner-line ones. Following the collection of experimental ρ and ρ * data, the inversion procedure requires iterating on R T , R 0 , and σ until a minimum-error solution to Eqs. 38 through 41 is found. Standard multivariate optimization algorithms may be used.
To compute the above shape factors for each iteration, V outer is fixed and I and V inner are computed for either the tangential or orthogonal case. To get I one can sum current flux over the entire region covered by boundary elements: is the area of boundary element i. The quantity in brackets will already be known from the calculations to determine influence coefficient A ii . To get V inner we likewise reuse terms already determined (see Eq. 36):

Numerical Results and Discussion
Input parameters.-Here we illustrate the use of the 2D and 3D models, giving one an idea of what happens for different values of input parameters. The four-line probe that is used has the following geometric parameters: L = 500 μm, p 1 = 4 μm, p 2 = 16 μm, p 3 = 24 μm, p 4 = 36 μm. The center-to-center distance between adjacent probe lines is therefore P = 20 μm. The baseline electrode sample is taken to have thickness T = 40 μm, current-collector dimensionless contact resistance R 0 = 1.5, and probe dimensionless contact resistance R T = 0.35. These values roughly correspond to what we have observed for commercial-quality LiCoO 2 cathodes 1 . These baseline values are used in each figure below unless specified differently (of course, for the 2D model L and R T are not used). The actual conductivity σ is not specified because the shape factors and normalized potentials are independent of σ.
One also must specify the size of the simulation domain (W x , W z ) and the number of Fourier terms (N , M) to use in order to maintain needed accuracy. Generally these limits need to increase with T and R 0 to account for the degree that current spreads out from the point of current injection. The models were run repeatedly and compared to cases where the domain size, number of Fourier terms, or both were set to large values. Based on that comparison, the following empirical relations were developed.
The floor function truncates the argument to the next smaller integer.
Recall that x = p 4 and z = L are the outer edges of the outer electrode, so the above relations are designed to place the domain boundaries W x and W z well past this point. These limits were found to result in relative errors in 2D shape factors (both tangential and orthogonal) of 0.03% or less for a wide range of thicknesses and contact resistances, namely 10 μm ≤ T ≤ 200 μm and 0 ≤ R 0 ≤ 4000. For corresponding 3D inner shape factors, relative errors were 0.1% or less for a similar wide range of parameters, including 100 μm ≤ L ≤ 1000 μm.
The two model solutions were implemented in C++ code. Results were also compared to 3D FEM solutions from COMSOL Multiphysics. For COMSOL to match the level of accuracy obtained by our 3D model required at least 10 5 tetrahedral elements, even with optimal sizing and distribution. With the range of input parameters we explored, a laptop computer (Intel i7 processor, 8GB RAM) could Figure 5. Fraction of current injected by the probe that travels through the current collector during the tangential experiment, as a function of sample thickness. Results for the 2D and 3D baseline models are virtually identical for this property. solve our 3D model (i.e. generate a shape factor) in 0.01 to 3 seconds. For the 2D model solution times were much less time than this. The comparable COMSOL FEM calculations each took 1 to 10 minutes. In addition, the FEM became problematic for some sets of input parameters due to extreme memory requirements in order to get sufficient accuracy. Thus, the FEM did not lead to a satisfactory and robust inversion technique.
Effect of current collector.-The main motivation for developing the four-line probe was to enable measurement of thin electrode films attached to the current collector. The model can predict how much shunt current during a tangential (i.e. traditional) experiment passes through the current collector. (For the orthogonal experiment by design all the current passes into the current collector.) Fig. 5 shows how as the thickness T of the sample decreases, the fraction of shunt current increases. The larger the fraction of shunt current, the more difficult it is to reliably determine σ from the tangential experiment, because the shunt current involves a modified path with resistances besides that of the bulk film. This result-and others below-could be generalized by thinking in terms of the ratio of T to the inter-line spacing P. This is the reason for making P small enough that most of the current bypasses the current collector. On the other hand, if P is too small then the probe may sample conductivity only near the surface rather than throughout the full cross section. In addition, particle-size effects can make the experiment less reliable if the contact widths and gaps are smaller than the largest particle sizes, which are around 10 μm for a typical Li-ion battery electrode.
Potential distributions. -Fig. 6 shows how potential varies at select locations inside the sample for the tangential and orthogonal experiments. Note that potentials are normalized for each plot. Parts a and b allow a comparison between typical tangential and orthogonal experiments. The maxima for the y = T curves occur, as expected, at the region of current injection, which is the outer line (24 μm ≤ x ≤ 36 μm). The vertical height of the dotted line (potential at y = 0) is proportional to the amount of local current density flowing into the current collector. As expected, in the case of the orthogonal experiment all of the current flows into the current collector. Fig. 6a allows one to compare the 2D model (R T = ∞) to the 3D model (R T = 0.35) for the tangential experiment. The change in the top boundary condition, for instance, results in modest perturbations to the potential in the vicinity of the inner line (4 μm ≤ x ≤ 16 μm). Though not shown in Fig. 6, the potential distribution for the 3D orthogonal experiment follows a similar trend. Fig. 6c shows the effect of a finite line length in the 3D model. The plot gives the potential in the sample running along the center of the inner contact line. The potential departs from the maximum as the end of the inner line is approached (z = L = 500 μm), indicating the spreading or leakage of current around the end of the contact lines (z > L). This spreading effect increases with T and R 0 values, as reflected in the need to increase W z under these conditions (Eq. 45). As given in Eq. 18, the potential is averaged to determine V inner . This suggests that V inner , and hence shape factors, for the 3D model will be less than those for the 2D model, to greater or lesser degree. Provided that R T is sufficiently large, the discrepancy between the 2D and 3D models is on the order of T /L. This is confirmed in Figs. 7 and 8. Fig. 7 shows how the inner shape factors change with electrode thickness and contact resistance. Figs. 7a and 7b present the tangential results for the 2D and 3D models, which are in moderate agreement. As T increases each model converges to a constant shape factor irrespective of R 0 . In the case of the 2D model this is given by an analytical result (Eq. 15), S t ∞ = 0.457 for the present geometry. Conversely, at low T values changing R 0 has a large effect on the tangential shape factor, because so much of the current is moving through the current collector (cf. Fig. 5) and would cause dependence on current collector contact resistance. Figs. 7c and 7d give corresponding shape factors for the orthogonal experiment. In this case there is a more notable difference between 2D and 3D results, meaning that the 2D result is less reliable. For the 2D model, shape factors increase without bound as T increases; for the 3D model, orthogonal shape factors asymptotically approach a constant value with increasing T , similar to the tangential case. The 3D asymptotes in Figs. 7b and 7d are complicated functions of the geometry, and no analytical value is known for the general case.

Shape factors.-
One compares S to experimental apparent resistivities ρ (Eqs. 38 through 41) to get the needed conductive properties of the sample. S and ρ values involve a voltage in the numerator and a current in the denominator. In typical electrical experiments, uncertainty in V can be high. If V inner is in the millivolt range, for instance, then relative error in ρ can be very high. In practical terms what this means is that any experiments conducted where S (and hence ρ) is too close to zero will be problematic and lead to large uncertainty in derived values such as σ. This will be apparent to an experimenter when an experiment is repeated to get information on the relative uncertainty in ρ. If S is too small to permit accurate results, it can be increased in many cases by decreasing the ratio of probe spacing P to sample thickness T . This means that the spacing between probe contacts must be carefully matched to the thickness of the sample-something which has not necessarily been considered in prior work. Fig. 8 shows how the inner shape factors S and outer shape factors S * change with R T , the resistance between the sample and probe. All parameters are kept at their baseline values except for R T . In addition, a second set of curves is generated with a thinner-than-baseline sample, T = 20 μm, to give an idea of the interplay between R T and thickness. As shown in Fig. 8a, the inner shape factors, particularly the orthogonal case, are relatively insensitive to R T over a wide range of values. We observe asymptotic behavior in S for very large or small R T values, with the crossover point around 0.1 to 1. Corresponding points (filled diamonds) are shown for the 2D model. In each case the 2D model overpredicts S by 2 to 5% relative to the 3D model. We do not expect perfect agreement, even when the 3D model uses a large R T value, as the 2D model also assumes L → ∞ and therefore neglects the leakage of current around the ends of the contact lines, as discussed above.
Theoretically if the areas of the four probe contacts were zero, then there would be no dependence of S on R T . Nevertheless, if one wants to generate the most accurate results, then R T must be accounted for by means of additional experimental information in the form of outer shape factors S * . These are given in Fig. 8b. In this case, S * is strongly dependent on R T when R T is large, in contrast to S in Fig. 8a. Therefore, the complementary combination of inner and outer line voltages (i.e. S and S * ) allows σ and R 0 -the parameters we care most about-to be determined accurately.

Conclusions
Measurement of electronic conduction properties of porous thinfilm battery electrodes poses significant challenges, particularly when the film is attached to a metallic current collector. We have developed a mathematical model of a four-line probe and a testing procedure that overcomes many of these difficulties while relying on principles similar to commonly used four-point probes. The procedure includes both tangential and orthogonal experiments to allow the determination of two properties: bulk electronic conductivity of the film and contact resistance with the current collector. Results are cast in terms of dimensionless shape factors that are functions of the geometry and the contact resistances. The ideal probe line spacing P is close to the value of the electrode film thickness T . Such an optimum involves tradeoffs between relative noise and bias in the voltage signal and the distribution of the volume sampled.
The 2D and 3D versions of the model were designed to allow rapid inversion of experimental data, in which shape factors must be computed multiple times in an iterative scheme. The 3D solution required a hybrid numerical/analytical method to obtain reasonable computational cost, though it is still more computationally demanding than the 2D solution. The new methods appear to be much more computationally efficient than corresponding FEM solutions, based on a modest set of input parameters that were tried. Moreover, there is scientific value in a more analytical approach. For example, determining shape factors at extremely large or extremely small resistance or conductivity values can cause problems in purely numerical solutions as the equations become stiff. Lastly, our models provide a general solution that can be tailored quickly to fit any probe or sample geometry. Indeed, the 3D model could easily be adapted to non-battery-related systems, or even a four-point probe in which the contact areas are arbitrary shapes.