Superconducting receiver arrays for magnetic resonance imaging

Sup erconducting receiv er arra ys for magnetic resonance imaging K o os C. J. Zev enho ven, ∗ An tti J. Mäkinen, and Risto J. Ilmoniemi Dep artment of Neur oscienc e and Biome dic al Engine ering, A alto University Scho ol of Scienc e, FI-00076 AAL TO, Finland Sup erconducting QUan tum-Interference Devices (SQUIDs) make magnetic resonance imaging (MRI) p ossible in ultra-lo w microtesla-range magnetic fields. In this w ork, we in vestigate the design parameters affecting the signal and noise p erformance of SQUID-based sensors and multic hannel magnetometers for MRI of the brain. Besides sensor in trinsics, v arious noise sources along with the size, geometry and n umber of sup erconducting detector coils are imp ortant factors affecting the image quality . W e derive figures of merit based on optimal combination of m ultichannel data, ana- lyze differen t sensor arra y designs, and pro vide tools for understanding the signal detection and the differen t noise mechanisms. The w ork forms a guide to making design decisions for b oth imaging- and sensor-oriented readers. CONTENTS I. In tro duction 1 I I. SQUID-detected MRI 2 A. Signal mo del and single-c hannel SNR 2 B. SQUIDs, pickup coils and detection 4 C. Sensitivit y patterns and signal scaling 4 I II. Noise mechanisms and scaling 6 A. Flux coupling and SQUID noise 6 B. Thermal magnetic noise from conductors 7 C. MRI electronics, coils and other noise sources 8 IV. Sensor arrays 9 A. Com bining data from multiple channels 9 B. Figures of merit and scaling for arrays 10 C. Correlation of noise b et ween sensors 10 D. Filling the array 11 V. Metho ds for numerical study 11 VI. Results 11 A. Effects of size and num ber 11 B. Realistic sensor configurations 12 VI I. Conclusions and outlo ok 13 References 13 App endix A. Exact Biot–Sav art integral ov er p olygonal path 16 I. INTR ODUCTION Magnetic r esonanc e imaging (MRI) is a widely used imaging metho d in clinical applications and research. It is based on measuring the magnetic signal resulting from ∗ koos.zevenho ven@aalto.fi nucle ar magnetic r esonanc e (NMR) of 1 1 H n uclei (pro- tons). In NMR, the magnetization rotates around an applied magnetic field ~ B at the proton Larmor frequency f L , whic h is prop ortional to B [1]. This b ehavior of the magnetization is often referred to as pr e c ession due to the direct connection to the quantum mechanical precession of nuclear spin angular momentum. Con ven tionally , the magnetic precession signal has b een detected using induction coils. The v oltage induced in a coil b y an oscillating magnetic field is prop ortional to the frequency of the oscillation, leading to v anishing signal amplitudes as f L approac hes zero. T o day , clin- ical MRI scanners indeed use a high main static field ~ B 0 ; t ypically B 0 = 3 T, corresponding to a frequency f 0 = 128 MHz. Ho wev er, when the signal is detected using magnetic field (or flux) sensors with a frequency- indep enden t resp onse, this need for high frequencies dis- app ears. Combined with the so-called prep olarization tec hnique for signal enhancement, highly sensitiv e mag- netic field detectors, t ypically those based on sup er c on- ducting quantum-interfer enc e devic es (SQUIDs), provide an NMR signal-to-noise ratio (SNR) that is indep endent of B 0 [2]. In recen t years, there has b een growing interest in ultra-low-field (ULF) MRI, usually measured in a field on the order of Earth’s magnetic field ( B 0 ∼ 10 – 100 µ T). A num b er of ULF-MRI-sp ecific imaging techniques ha ve emerged, including rotary-scanning acquisition (RSA) [4], temp erature mapping [5], signal-enhancing dynamic nuclear p olarization [6, 7], imaging of electric curren t density (CDI) [8–10], and making use of sig- nifican t differences in NMR relaxation mec hanisms at ULF compared to tesla-range fields [11–13]. Sev eral groups hav e also in vestigated p ossibilities to directly de- tect c hanges in the NMR signal due to neural currents in the brain [14–17] and electrical activ ation of the heart [18]. A further notable field of research now fo cuses on combining ULF MRI with magneto encephalography (MEG). In MEG, an arra y of typically ∼ 100 sensors [19– 21] is arranged in a helmet-shap ed configuration around the head (see Fig. 1) to measure the weak magnetic fields pro duced by electrical activity in the brain [22, 23]. SQUID sensors tailored for ULF MRI can typically also 2 Figure 1. Helmet-type sensor array geometries consisting of (a) triple-sensor mo dules at 102 p ositions similar to stan- dard Elekta/Neuromag MEG configurations and (b) an arra y with larger ov erlapping pickup coils for increased p erfomance. Magnetometers are marked in green and gradiometers in red or blue; see Sec. I I B for descriptions of pickup coils. (The sample head shap e is from MNE-Python [3].) b e used for MEG, and p erforming MEG and MRI with the same device can significantly impro ve the precision of lo calizing brain activit y [24–28]. In typical early ULF-MRI setups [2], the signal w as detected b y a single dc SQUID coupled to a sup ercon- ducting pickup coil w ound in a gradiometric configura- tion that rejects noise from distant sources. In this case, the maximum size of the imaging field of view (F OV) is roughly given b y the diameter of the pickup coil. With large diameters such as 60 mm, field sensitivities b et- ter than 1 fT / √ Hz hav e b een ac hieved with a reason- able F OV. A large coil size, how ever, do es ha ve its dra w- bac ks, including issues such as high inductance and in- creased requirements in dynamic range. Therefore, the most straightforw ard wa y to increase the av ailable FO V and the SNR is to use an array of sensors. In addition, as is well kno wn in the context of MEG [20, 29, 30], a m ulti-channel measuremen t allows forming so-called soft- w are gradiometers and more adv anced signal pro cessing tec hniques to reduce noise that can be optimized sep- arately for differen t noise environmen ts. In ULF-MRI, this can even b e done individually for eac h v oxel (vol- ume element) p osition within the imaging target, as will b e sho wn later. While single-c hannel systems are still common, several groups hav e already b een using arra ys of sensors. Also in con v entional MRI, so-called parallel MRI is p erformed using an array of tens of induction coils, al- lo wing full reconstruction of images from a reduced num- b er of data acquisitions [31, 32]. There are studies on designing arrays of induction coils for parallel MRI [33] with an emphasis on minimizing artefacts caused by the reduced n umber of acquisitions. At the kHz frequencies of ULF MRI, the dominant noise mechanisms are signif- ican tly different, and one needs to consider, for instance, electromagnetic interference from p o wer lines and electri- cal equipmen t, thermal noise from the radiation shield of the cryostat required for op erating the sup erconducting sensors, as well as noise and transients from other parts of the ULF MRI system structure and electronics [34]. Studies on the design of arrays for MEG [20, 35, 36], whic h mainly fo cus on the accuracy of lo calizing brain activit y , are also not applicable to ULF MRI. In terms of single-sensor ULF-MRI signals, there are existing stud- ies of the depth sensitivity [37] and SNR as a function of frequency with different detector types [38]. Previously , in Ref. [39], we presen ted approaches for quan titative comparison of sensor arrays in terms of the com bined performance of the sensors, the results indicat- ing that the optimum sensor for ULF MRI of the brain w ould b e somewhat larger than t ypical MEG sensors. Extending and refining those studies, we aim to provide a fairly general study of the optimization of ULF-MRI arra y p erformance, with sp ecial attention to SNR and imaging the human head. W e b egin by defining relev ant quantities and review- ing basic principles of ULF MRI in Sec. I I. Then, w e analyze the effects of sensor geometry and size with dif- feren t noise mechanisms (Sec. I I I), adv ancing to sensor arra ys (Sec. IV). Finally , we show computed estimations of array SNR as functions of pic kup size and num b er, and pro vide more detailed comparison of spatial SNR profiles with different array designs (Secs. V and VI). I I. SQUID-DETECTED MRI A. Signal mo del and single-channel SNR In contrast to conv entional MRI, where the tesla-range main field is static and accounts for b oth polarizing the sample and for the main readout field, ULF MRI employs switc hable fields. Dedicated electronics [34] are able to ramp on and off even the main field ~ B 0 with an ultra-high effectiv e dynamic range. An additional pulsed prep olar- izing field ~ B p magnetizes the target b efore signal acqui- sition. Typically , a dedicated coil is used to generate ~ B p ( B p ∼ 10 – 100 mT) in some direction to cause the proton bulk magnetization ~ M ( ~ r ) to relax with a longitudinal re- laxation time constant T 1 to wards its equilibrium v alue corresp onding to ~ B p . After a p olarizing time on the or- der of seconds or less, ~ B p is switc hed off—adiabatically , in terms of spin dynamics—so that ~ M turns to the direc- tion of the remaining magnetic field, t ypically ~ B 0 , while k eeping most of its magnitude. Next, say at time t = 0 , a short excitation pulse ~ B 1 is applied whic h flips ~ M a wa y from ~ B 0 , typically by 90 ◦ , bringing ~ M in to precession around the magnetic field at p ositions ~ r throughout the sample. While rotating, ~ M ( ~ r ) decays to wards its equilibrium v alue corresp onding to the applied magnetic field in which the magnetization precesses. This field, ~ B L , ma y sometimes simply b e a uniform ~ B 0 , but for spatial enco ding and other purp oses, differen t non-uniform magnetic fields ∆ ~ B ( ~ r, t ) are addi- tionally applied to affect the precession before or during 3 acquisitions. The enco ding is tak en in to accoun t in the subsequen t image reconstruction. The ULF MRI signal can b e modeled to a high accu- racy giv en the absence of unstable distortions common at high frequencies and high field strengths. T o obtain a mo del for image formation, we b egin b y examining ~ M at a single p oint. If the z axis is set parallel to the total precession field ~ B L , then the xy (transverse) comp onen ts of ~ M accoun t for the precession. Assuming, for no w, a static ~ B L , and omitting the decay for simplicity , the transv erse magnetization ~ M xy = ~ M xy ( t ) can be written as ~ M xy ( t ) = M xy [ b e x cos( ω t + φ 0 ) − b e y sin( ω t + φ 0 )] , (1) where ω = 2 π f L is the precession angular frequency , b e ♥ is the unit v ector along the ♥ axis ( ♥ = x, y , z ), and φ 0 is the initial phase, which sometimes contains useful information. In an infinitesimal v olume dV at p osition ~ r in the sam- ple, the magnetic dip ole momen t of protons in the v ol- ume is ~ M ( ~ r ) dV . It is straightforw ard to show that the rotating comp onen ts of this magnetic dip ole are seen by an y magnetic field or flux sensor as a sinusoidal signal dψ s = | β | cos( ω t + φ 0 + φ s ) M xy dV . Here | β | = | β ( ~ r ) | is the p eak sensitivity of the sensor to a unit dip ole at ~ r that precesses in the xy plane, and φ s = φ s ( ~ r ) is a phase shift depending on the relative p ositioning of the sensor and the dip ole. T o obtain the total sensor signal ψ s , dψ s is integrated ov er all space: ψ s ( t ) = Z | β ( ~ r ) | M xy ( ~ r ) cos φ ( ~ r , t ) d 3 ~ r , (2) where φ ( ~ r, t ) = Z t 0 ω ( ~ r , t 0 ) dt 0 + φ 0 ( ~ r ) + φ s ( ~ r ) . Here, w e hav e noted that the magnetic field can v ary in b oth space and time and therefore ω = ω ( ~ r , t ) = γ B ( ~ r, t ) , where γ is the gyromagnetic ratio; γ / 2 π = 42 . 58 MHz/T for a proton. F or con venience, the signal given by Eq. (2) can be demo dulated at the angular Larmor frequency ω 0 = 2 π f 0 corresp onding to B 0 ; using the quadrature comp onent of the phase sensitive detection as the imaginary part, one obtains a complex-v alued signal Ψ ( t ) = Z | β ( ~ r ) | M xy ( ~ r ) e − i [ φ ( ~ r,t ) − ω 0 t ] d 3 ~ r = Z β ∗ ( ~ r ) m ( ~ r ) e − i R t 0 ∆ ω ( ~ r ,t 0 ) dt 0 d 3 ~ r , (3) where ∗ denotes the complex conjugate, m ( ~ r ) = M xy ( ~ r ) e − iφ 0 ( ~ r ) is the uniform-sensitivity image , ∆ ω = ω − ω 0 , and we define β ( ~ r ) = | β ( ~ r ) | e iφ s ( ~ r ) (4) as the single-c hannel c omplex sensitivity pr ofile . Besides geometry , β generally also dep ends on the direction of the precession field; β = β ~ B L ( ~ r ) . After acquiring enough data of the form of Eq. (3), the image can b e reconstructed—in the simplest case us- ing only one sensor, or using multiple sensors, each ha v- ing its o wn sensitivity profile β . As a simplified mo del for understanding image formation, ideal F ourier enco d- ing turns Eq. (3) into the 3-D F ourier transform of the sensitivit y-weigh ted complex image β ∗ m = ( β ∗ m )( ~ r ) . In realit y , how ev er, the inv erse F ourier transform only pro- vides an approximate reconstruction, and more sophisti- cated techniques should b e used instead [40]. Here, w e do not assume a sp ecific spatial enco ding sc heme. Notably , how ever, the sensitivity profile is in- distinguishable from m based on the signal [Eq. (3)]. In other w ords, the spatial v ariation of β ∗ affects the ac- quired data in the same wa y as a similar v ariation of the actual image w ould, regardless of the spatial enco ding sequence in ∆ ω . Consider a small vo xel of centered at ~ r . The contribu- tion of the vo xel to the signal in Eq. (3) is prop ortional to an effective v oxel v olume V . Due to measurement noise, the vo xel v alue b ecomes V β ∗ m + ξ , where ξ is a random complex noise term. If β is known, the in tensit y- corrected v oxel of a real-v alued image from a single sensor is given by Re  m ( ~ r ) + ξ V β ∗ ( ~ r )  = m ( ~ r ) + Re  ξ e iφ s  | s ( ~ r ) | , (5) where s ( ~ r ) = V β ∗ ( ~ r ) is the sensitivit y of the sensor to m in the giv en vo xel. Assuming that the distribution of ξ = | ξ | e iφ ξ is indep endent of the phase φ ξ , the standard deviation σ of Re  ξ e iφ s  is indep enden t of φ s and pro- p ortional to σ s , the standard deviation of the noise in the relev ant frequency band of the original sensor signal. The precision of a vo xel v alue can b e describ ed by the (amplitude) SNR of the vo xel v alue. The vo xel SNR is defined as the correct vo xel v alue m ( ~ r ) divided by the standard deviation of the random error and can b e writ- ten as SNR = m ( ~ r ) V | β ( ~ r ) | σ ∝ B p V | β ( ~ r ) | √ T tot σ s , (6) where the last expression incorp orates that m ∝ B p , and that σ is inv ersely prop ortional to the square ro ot of the total signal acquisition time, whic h is proportional to the total MRI scanning time T tot . It should b e recognized, ho wev er, that σ also dep ends heavily on factors not visi- ble in Eq. (6), such as the imaging sequence. Ultimately , the abilit y to distinguish b et w een differ- en t types of tissue depends on the c ontr ast-to-noise r atio (CNR), which can b e defined as the SNR of the differ- ence b et w een image v alues corresp onding to tw o tissues. A b etter CNR can b e ac hieved b y improving either the SNR or the contrast, which b oth strongly dep end also on the imaging sequence. 4 B. SQUIDs, pickup coils and detection SQUIDs are based on sup er c onductivity , the phe- nomenon where the electrical resistivity of a material completely v anishes b elow a critical temp erature T c [41]. A commonly used material is niobium (Nb), which has T c = 9 . 2 K. It is usually co oled by immersion in a liquid helium bath that b oils at 4 . 2 K in atmospheric pressure. SQUIDs can b e divided in to t w o categories, rf and dc SQUIDs, of whic h the latter is t ypically used for bio- magnetic signals as w ell as for ULF MRI [19, 27]. The dc SQUID is a sup erconducting lo op interrupted by tw o w eak links, or Josephson junctions; see Fig. 2(a). With suitable shun ting and biasing to set the electrical operat- ing p oint, the current or v oltage across the SQUID can b e configured to exhibit an oscillatory dep endence on the magnetic flux going through the loop—analogously to the w ell known double-slit interference of wa ves. A linear resp onse to magnetic flux is obtained by op- erating the SQUID in a flux-lock ed lo op (FLL), where an electronic control circuit aims to keep the flux constant b y applying negative flux feedbac k via an additional feed- bac k coil. T o av oid harmful resonances and to achiev e low noise, the SQUID lo op itself is usually made small. The sig- nal is coupled to it using a larger pickup coil connected to the SQUID via an input circuit to achiev e high sen- sitivit y . An input circuit may simply consist of a pickup c oil and an input c oil in series, forming a con tin uous su- p erconducting path whic h, by ph ysical nature, conserves the flux through itself, and feeds the SQUID according to the signal received by the pickup coil, as explained in Sec. I II A along with more sophisticated input circuits. Differen t types of resp onses to magnetic fields can b e ac hieved b y v arying the pickup coil geometry . Fig. 2(b–g) sc hematically depicts some p opular t yp es. The simplest case is just a single loop, a magnetometer , which in a ho- mogeneous field responds linearly to the field component p erpendicular to the plane of the lo op (b). T wo lo ops of the same size and orientation, but w ound in opp osite directions, can b e used to form a gr adiometer . The re- sulting signal is that of one lo op subtracted from that of the other. It can b e used to appro ximate a deriv ative of the field comp onen t with resp ect to the direction in whic h the lo ops are displaced (by distance b , called the baseline). Typical examples are the planar gradiome- ter (c) and the axial gradiometer (d). By using more lo ops, one can measure higher-order deriv atives. Some ULF-MRI implementations [2, 42] us e second-order axial gradiometers (e). If a source is close to one lo op of a long- baseline gradiometer, that ‘pickup lo op’ can b e thought of as a magnetometer, while the additional lo ops sup- press noise from MRI coils or distant sources. Ho wev er, adding lo ops also increases the inductance L p . Before a more detailed theoretical discussion regarding L p and SQUID noise scaling, we study the detection of the MRI signal by the pickup coils. Figure 2. Sc hematic (a) of a simple SQUID sensor and the flux-lo c k ed loop (more detail in Secs. I II A and IV C), and (b–f ) of differen t types of pickup coils. Pic kup coil types are (b) magnetometer (M0), (c) planar first-order gradiome- ter (PG1), (d) axial first-order gradiometer (AG1), (e) axial second-order gradiometer (AG2), (f ) planar gradiometer with a long baseline, and (g) a magnetometer and tw o planar gra- diometers in a triple-sensor unit (M0, PG1 x , PG1 y ). C. Sensitivity patterns and signal scaling The magnetic flux Φ pick ed up by a coil made of a thin sup erconductor is given by the integral of the magnetic 5 field ~ B ov er a surface S b ound by the coil path ∂ S , Φ = Z S ~ B · d 2 n ~ r = I ∂ S ~ A · d~ r . (7) Here, the line integral form was obtained b y writing ~ B in terms of the vector p otential ~ A as ~ B = ∇ × ~ A , and applying Stokes’s theorem. As explained in Sec. I I A, the signal in MRI arises from spinning magnetic dip oles. The quasi-static ap- pro ximation holds well at signal frequencies, providing a vector p oten tial for a dip ole ~ m positioned at ~ r 0 as ~ A ( ~ r ) = µ 4 π ~ m × ( ~ r − ~ r 0 ) | ~ r − ~ r 0 | 3 , where µ is the p ermeability of the medium, assumed to b e that of v acuum; µ = µ 0 . Sub- stituting this into Eq. (7) and rearranging the resulting scalar triple pro duct leads to Φ = ~ m · ~ B s ( ~ r 0 ) , ~ B s ( ~ r 0 ) = µ 4 π I ∂ S d ~ r × ( ~ r 0 − ~ r ) | ~ r 0 − ~ r | 3 , (8) where the expression for the sensor field ~ B s is the Biot– Sa v art formula for the magnetic field at ~ r 0 caused by a h yp othetical unit current in the pickup coil, as required b y recipro city . The sensor field ~ B s is closely related to the complex sensitivit y pattern β introduced in Sec. I I A. In an applied field ~ B L = B L b e z , the magnetization precesses in the xy plane, and β can in fact b e written as β ( ~ r ) = ~ B s ( ~ r ) · ( b e x + i b e y ) . (9) F or arbitrary ~ B = B L b e L , we hav e | β ~ B ( ~ r ) | = q | ~ B s ( ~ r ) | 2 − [ ~ B s ( ~ r ) · b e L ] 2 . (10) W e choose to define the measured signal as the flux through the pic kup coil—a con ven tion that app ears throughout this pap er. The measuremen t noise is consid- ered accordingly , as flux noise. This contrasts looking at magnetic-field signals and noise, as is often seen in the lit- erature. W orking with magnetic flux signals allows for di- rect comparison of differen t pickup coil types. Moreo ver, the appro ximation that magnetometer and gradiometer pic kups resp ond to the field and its deriv atives, resp ec- tiv ely , is not alwa ys v alid. The signal often scales as simple p o wer laws R α with the pickup coil size R (or radius, for circular coils). When the distance l from the coil to the signal source is large compared to R , a magnetometer sees a flux Φ ∝ B R 2 , giving an amplitude sc aling exp onent α = 2 . When scal- ing a gradiometer, ho wev er, also the baseline b is pro- p ortional to R . This leads to α = 3 for a first-order gradiometer, or α = 2 + k for one of k th order. Con- v ersely , the signal scales with the distance as l − α − 1 , as is v erified by writing the explicit forms of the field and its deriv ativ es. The additional − 1 in the exp onen t re- flects the dip olar nature of the measured field ( − 2 for quadrup oles etc.). Figure 3. Isosurfaces of sensitivity patterns | β ( ~ r ) | inside a helmet array for tw o of the magnetometer lo ops mark ed in red. The arrow depicts the direction of the precession field ~ B L during readout ( e.g. ~ B 0 ). Note that, b ecause of the precession plane, there are insensitiv e directions (“blind angles”) in the profiles, dep ending on the relative orien tation of ~ B L . F or some cases, the detected flux can be calculated analytically using Eq. (8). First, as a simple example, consider a dipole at the origin, and a circular magne- tometer pickup lo op of radius R parallel to the xy plane at z = l , centered on the z axis. The in tegral in Eq. (8) is easily integrated in cylindrical co ordinates to give ~ B s = B s b e z = µR 2 2( R 2 + l 2 ) 3 2 b e z . (11) If the dip ole precesses in, for instance, the xz plane, the corresp onding sensitivity is | β | = B s . Instead, if preces- sion takes place in the xy plane, the sensitivity v anishes; | β | = 0 , and no signal is receiv ed. In this case, moving the pickup lo op aw a y from the z axis w ould cause a sig- nal to app ear. These extreme cases show that even the absolute v alue of a single-c hannel sensitivit y is strongly dep enden t on the sensor orien tation with respect to the source and the magnetic field, as is also seen in Fig. 3. Another notable prop erty of the sensitivity | β | = B s from Eq. (11) is that if l is fixed, there is a v alue of R 6 ab o v e which the sensitivity starts to decrease, i.e. , part of the flux going through the lo op comes back at the edges canceling a portion of the signal. By requiring ∂ B s /∂ R to v anish, one obtains R = l √ 2 , the lo op radius that gives the maximum signal. In terestingly , how ever, if instead of the perp endicular ( z ) distance, l is taken as the closest distance to the pickup-coil winding, then the coil is on a spherical surface of radius R a = l . No w, based on Pythagoras’s theorem, R 2 + l 2 in Eq. (11) is replaced with l 2 . In other words, the sensor field is simply ~ B s = b e z µR 2 / 2 l 3 , so the scaling of α = 2 happ ens to b e the same as for distant sources in this simple case. Imp ortan tly , how ever, the noise mec hanisms also de- p end on R , and moreov er, the situation is complicated b y the presence of multiple sensors. These matters are discussed in Secs. I I I–IV. I II. NOISE MECHANISMS AND SCALING The signal from each measurement channel, corre- sp onding to a pickup coil in the sensor arra y , contains flux noise that can originate from v arious sources. Examples of noise sources are the sensor itself, noise in electron- ics that drives MRI coils, cryostat noise, magnetic noise due to thermal motion of particles in other parts of the measuremen t device and in the sample, noise from other sensors, as well as environmen tal noise. This section is dev oted to examining the v arious noise mechanisms and ho w the noise can b e dealt with. Unless stated otherwise, noise is considered a random signal with zero a verage. W e use amplitude scaling exp onen ts α to characterize the dep endence of noise on pickup-coil size and type. A. Flux coupling and SQUID noise F or estimates of SQUID sensor noise as a function of pic kup coil size, a mo del for the sensor is needed. As explained in Sec. I I B, the signal is coupled into the SQUID loop via an input circuit. In general, the in- put circuit may consist of a sequence of one or more all- sup erconductor closed circuits connected by intermediate transformers. Via inductance matc hing and coupling op- timization, these circuits are designed to efficien tly couple the flux signal into the SQUID lo op. In termediate transformers can be useful for optimal coupling of a large pic kup coil to a SQUID-coupled input coil, as analyzed e.g. in Ref. [43]. T o further understand the concept, consider a tw o-stage input circuit where a pic kup coil ( L p ) is connected to a transmitting inductor L 1 to form a closed superconducting path; see Fig. 4. Ideally , the distance b etw een the tw o coils is fairly small in order to av oid signal loss due to parasitic inductances of the connecting traces or wiring. The total inductance of this flux-coupling circuit by itself is L p + L 1 . The pri- mary is coupled to a secondary inductor L 2 with mutual inductance M 12 . As the magnetic flux pick ed up in L p Figure 4. Simplified schematic of a superconducting SQUID input circuit. Zero or more in termediate transformers (dashed b o x) may be present. c hanges by ∆ Φ p , there is a corresp onding change ∆ J 1 in the sup ercurren t flowing in the circuit such that the flux through the closed path remains constant. This passes the flux signal onw ards to L 2 whic h forms another flux- transfer circuit together with the input coil L i , which couples inductively into the SQUID. Sup erconductivit y has tw o imp ortant effects on the transmission of flux in to the next circuit. First, the pres- ence of sup erconducting material close to a coil tends to reduce the coil inductance because of the Meissner effect: the magnetic flux is exp elled and the material acts as a p erfect diamagnet. This effect is included in the given in- ductances L p and L 1 . The other effect emerges when the flux is transmitted in to another closed sup erconducting circuit, such as via M 12 . This is b ecause the transmitting coil is sub ject to the counteracting flux M 2 12 ∆ J 1 / ( L 2 + L i ) from the receiving coil of the other circuit. No w curren t ∆ J 1 only generates a flux [ L 1 − M 2 12 / ( L 2 + L i )]∆ J 1 in L 1 . Closing the secondary circuit th us changes the in- ductance from L 1 to L 0 1 = L 1 − M 2 12 L 2 + L i = L 1  1 − k 2 12 1 + L i /L 2  , (12) where the last form is obtained by expressing the m u- tual inductance in terms of the coupling constant k 12 ( | k 12 | < 1 ) as M 12 = k 12 √ L 1 L 2 . Note that we do not include a coun teracting flux from the SQUID inductance L S bac k into L i , i.e. , no screening from the biased SQUID lo op. How ever, like other inductances, L i do es include the effect of the presence of the nearb y sup erconductors through the Meissner effect. The change of flux though the dc SQUID lo op is now obtained as ∆ Φ S = M iS ∆ J 2 = M iS M 12 L 2 + L i ∆ J 1 (13) = M iS M 12 ( L 2 + L i )( L p + L 1 ) − M 2 12 ∆ Φ p , (14) or, with M iS = k iS √ L i L S and defining χ 1 and χ 2 suc h that L 1 = χ 1 L p and L 2 = χ 2 L i , we hav e ∆ Φ S ∆ Φ p = k iS √ L S p L p × k 12 √ χ 1 χ 2 χ 1 χ 2 (1 − k 2 12 ) + χ 1 + χ 2 + 1 . (15) 7 F or a given pickup coil, χ 1 and χ 2 can usually b e cho- sen to maximize the flux seen b y the SQUID. While the function in Eq. (15) is monotonous in k 12 , there is a single maxim um with resp ect to parameters χ 1 , χ 2 > 0 . Noting the symmetry , we must hav e χ 1 = χ 2 =: χ , and the fac- tor in Eq. (15) b ecomes k 12 χ/ [ χ 2 (1 − k 2 12 ) + 2 χ + 1] , which is maximized at χ = 1 / p 1 − k 2 12 . At the optimum, the coupled flux is given by ∆ Φ S ∆ Φ p = k iS k 12 √ L S 2 p L p  1 + p 1 − k 2 12  − → k 12 → 1 − k iS 2 s L S L p . (16) Notably , with a k 12 ≈ 1 , the coupling corresp onds to a p erfectly matched single flux-coupling circuit [41]. Al- ready at k 12 = 0 . 8 , 50% of the theoretical maxim um is ac hieved, while matc hing without an intermediate trans- former may cause practical difficulties or parasitic reso- nances. When referred to SQUID flux Φ S , the noise in the mea- sured SQUID voltage in the flux-lo ck ed lo op corresp onds to a noise sp ectral density S Φ S ( f ) at frequency f . As the signal transfer from the pic kup coil to the SQUID is giv en b y Eqs. (15), the equiv alen t flux resolution referred to the signal through the pickup coil can be written as S 1 / 2 Φ p ( f ) = 2 p L p  1 + p 1 − k 2 12  k iS k 12 √ L S S 1 / 2 Φ S ( f ) . (17) Due to resonance effects and thermal flux jumps, L S needs to b e k ept small [41]. The flexibility of intermedi- ate transformers allo ws the same model to estimate noise lev els with a wide range of pickup coil inductances L p . In general, the inductance of a coil with a giv en shap e scales as the linear dimensions, or radius R , of the coil. If the wire thic kness is not scaled accordingly , there will b e an extra logarithmic term [44]. Even then, within a range small enough, the dep endence is roughly S 1 / 2 Φ p ∝ R α with α = 1 / 2 . The case of a magnetometer lo op in a homogeneous field then still has a field resolution S 1 / 2 B ( f ) prop ortional to R − 3 / 2 . B. Thermal magnetic noise from conductors Electric noise due to the thermal motion of charge car- riers in a conducting medium is called Johnson–Nyquist noise [45, 46]. According to Amp ` e re’s law ∇ × ~ B = µ 0 ~ J , the noise currents in the curren t densit y ~ J also produce a magnetic field which may in terfere with the measure- men t. In this view, devices should b e designed in such a w ay that the amount of conducting materials in the vicinit y of the sensors is small. How ever, there is a low er limit set by the conducting sample—the head. Estima- tions of the sample noise [38] hav e giv en noise lev els b e- lo w 0 . 1 fT / √ Hz , consisten t with a recent exp erimen tal result of 55 aT / √ Hz [47]. Other noise sources still ex- ceed those v alues b y more than an order of magnitude. More restrictingly , it is difficult to av oid metals in most applications. T o keep the SQUID sensors in the sup erconducting state, the array is k ept in a helmet-b ottom cryostat filled with liquid helium at 4 . 2 K. The thermal sup erinsulation of a cry ostat usually in volv es a v acuum as well as lay ers of aluminized film to suppress heat transfer by radia- tion [41]. The magnetic noise from the sup erinsulation can b e reduced by breaking the conducting materials into small isolated patches. Seton et al. [48] used aluminium- coated p olyester textile, whic h efficien tly breaks up cur- ren t paths in all directions. By using very small patc hes, one can decrease the field noise at the sensors b y orders of magnitude, although with increased He b oil-off [49]. T o lo ok at the thermal noise from the insulation lay ers in some more detail, consider first a thin slab with con- ductivit y σ on the xy plane at temp erature T . Johnson– Nyquist currents in the conductor pro duce a magnetic field ~ B ( x, y , z , t ) outside the film. F or an infinite (large) slab, the magnitude of the resulting field noise dep ends, b esides the frequency , only on z , the distance from the slab (assume z > 0 ). A t low frequencies, the sp ectral densities S B α ( α = x, y , z ) corresp onding to Cartesian field noise comp onen ts are then given by [50] S 1 / 2 B z = √ 2 S 1 / 2 B x = √ 2 S 1 / 2 B y = µ 2 s k B T 2 π σ d z ( z + d ) , (18) where d is the thickness of the slab and k B the Boltzmann constan t. The infinite slab is a goo d approximation when using a flat-bottom cryostat or when the radius of curv ature of the cryostat wall is large compared to individual pic kup lo ops. Consider a magnetometer pickup lo op with area A placed parallel to the conducting films in the insulation— to measure the z component of the magnetic field, B z . The coupled noise flux is the integral of B z o ver the lo op area. If the lo op is small, the noise couples to the pickup circuit as S 1 / 2 Φ = S 1 / 2 B z A . A coil of size R then sees a flux noise prop ortional to S 1 / 2 B z R 2 , that is, α = 2 . Instead, if the pickup coil is large, the situation is quite differen t. The instantaneous magnetic field dep ends on all co ordinates and v aries significantly o ver the large coil area. Consider the noise field at tw o points in the plane of the coil. The fields at the tw o p oints are nearly equal if the p oin ts are close to eac h other. Ho w ever, if the p oin ts are separated by a distance larger than a corre- lation length λ c ( z ) , the fields are uncorrelated. There- fore, if R  λ c , the coupled flux is roughly a sum of A/λ 2 c uncorrelated terms from regions in whic h the field is correlated. Eac h term has a standard deviation of or- der S 1 / 2 B z λ 2 c . The sp ectral density of the cry ostat noise is then S Φ, c ( f ) ≈ AS B z ( ~ r , f ) λ 2 c ( ~ r ) . (19) Most imp ortan tly , the flux noise amplitude S 1 / 2 Φ, c is di- rectly prop ortional to the coil size R , and w e now hav e 8 α = 1 . Still, the noise increases to a higher p o wer of R than the sensor noise, which according to section I II A scales as √ R and hence dominates in small pickup coils. F or a contin uous film, the correlation length λ c can b e estimated from data in Ref. [51] to b e around sev eral times z . The correlation at distances smaller than λ c is due to tw o reasons. First, the magnetic field due to a small current element in the conductor is spread in space according to the Biot–Sa v art law. Second, the noise cur- ren ts in elemen ts close to each other are themselves cor- related. The latter effect is broken down when the film is divided into small patc hes; only v ery small current lo ops can o ccur, and the noise field starts to resem ble that of Gaussian uncorrelated magnetic p oin t dip oles through- out the surface. In this case, Eq. (18) is no longer v alid, but the appro ximate relation of Eq. (19) still holds—now with a smaller λ c . The magnetometer case is easily extended to first-order planar gradiometers parallel to the sup erinsulation lay ers [Fig. 2( b , f )]. F or a very small baseline, b  λ c , the field noise is effectively homogeneous and th us cancels out. Ho wev er, when b  λ c , the spectral densit y of the noise p o w er is twice that of a single lo op. C. MRI electronics, coils and other noise sources As explained in Sec. I I A, MRI mak es heavy use of applied magnetic fields. The fields are generated with dedicated curren t sources, or amplifiers, to feed currents in to coils wound in different geometries. As opposed to applying static fields, a ma jor challenge arises from the need for oscillating pulses and the desire to quickly switc h on and off all fields, including not only readout gradien ts but also the main field ~ B 0 , whic h requires an ultra-high dynamic range to a void excess noise. Switching of ~ B 0 en- ables full 3-D field mapping for imaging of small electric curren ts in volume [34]. Noise in the coil curren ts can b e a ma jor concern in the instrumen tation. The con- tribution from ~ B 0 ideally scales with pickup coil size as R α , α = 2 for a magnetometer, and noise in linear gradi- en ts essentially scales as α = 2 in magnetometers as well as fixed-baseline gradiometers. With b ∝ R , first-order gradiometers exp erience noise from linear gradient coils according to α = 3 . MRI coils themselv es also pro duce Johnson–Nyquist noise. In particular, the p olarizing coil is often close to the sensors and made of thick wires as it should b e able to pro duce relativ ely high fields. This allo ws ther- mal electrons to form current lo ops that generate field noise with complicated spatial c haracteristics, whic h is detrimen tal to image quality and should b e eliminated. Another approach is to use litz wire, which is composed of thin wires individually coated with an insulating la yer. This prev ents significant noise currents p erp endicular to the wire and eliminates large current lo ops. Ho wev er, ef- ficien t uniform cooling of litz wire is problematic, leading to larger coil diameters. Increasing the coil size, how ever, significan tly increases harmful transients in the system as w ell as the pow er and co oling requirements [52]. Instead, w e hav e had promising results with thin custom-made sup erconducting filament wire and DynaCAN (Dynami- cal Coupling for Additional dimeNsions) in-sequence de- gaussing wa v eforms to solve the problem of trapp ed flux [52, 53]; optimized oscillations at the end of a pulse can exp el the flux from the sup erconductor. Such coils con- tain muc h less metal, and significantly reduce the size of curren t lo ops that can generate magnetic noise. A significant amount of noise also originates from more distan t lo cations. P ow er lines and electric devices, for instance, are sources that often can not b e remov ed. Indeed, magnetically shielded ro oms (MSRs) effectiv ely atten uate such magnetic interference. Ho wev er, pulsed magnetic fields inside the shielded ro om induce eddy curren ts exceeding 1 kA in conductiv e MSR walls [54], leading to strong magnetic field transients that not only saturate the SQUID readout, but also seriously interfere with the nuclear spin dynamics in the imaging field of view. Even a serious eddy current problem can again b e solv ed with a DynaCAN approach where optimized cur- ren t w av eforms are applied in additional coil windings to couple to the complexity of the transient [55]. Noise from distant sources t ypically scales with the pic kup coil size with an exp onen t at least as large as the signal from far-a wa y sources: α = 2 + k for a k th -order gradiometer (see Sec. II B). Although the noise detected b y gradiometers scales to a higher p o wer than with mag- netometers ( k = 0 ), gradiometers ha ve the adv antage that they , in principle, do not resp ond to a uniform fie ld. F or a higher-order gradiometer that is not to o large, the en vironmental noise is nearly uniform in space, and there- fore effectively suppressed by the pickup coil geometry . Gradiometers with relatively long baselines can also b e seen as magnetometers when the source is close to one of the lo ops. Still, they function as gradiometers from the p erspective of distant noise sources. A similar result ap- plies for so-called soft ware gradiometers, which can, for example, b e formed b y afterwards taking the difference of the signals of tw o parallel magnetometers. How ever, in Sec. IV A, a more sophisticated technique is describ ed for minimizing noise in the combination of multiple c hannels. A t very lo w system noise levels, other significant noise mec hanisms include noise due to dielectric los ses. Elec- trical activit y in the brain can also b e seen as a source of noise. This noise, ho wev er, is strongest at frequencies w ell b elo w 1 kHz. Using Larmor frequencies in the kHz range m a y therefore b e sufficien t for sp ectral separation of brain noise from MRI. The amplitude scaling exp onen ts α for signal and noise are summarized in T able I. The notation in later sections refers to the scaling of flux signal and noise in terms of α s and α n , respectively . F or a single sensor, the SNR scaling R δ is given by δ = α s − α n . 9 T able I. Amplitude scaling exp onen ts α for the flux noise stan- dard deviation σ ∝ R α as w ell as the signal, given different pic kup-coil geometries and noise mechanisms. Pic kup type (see Fig. 1) → M0 AG k PG k Sensor noise (optimally matched) 1/2 1/2 1/2 Sensor noise (unmatched, large L p ) 1 1 1 Distan t source, b ∝ R 2 2 + k 2 + k Distan t source, b fixed 2 2 – ~ B 0 amplifier 2 0 ∗ 0 ∗ Gradien t amplifiers, b ∝ R , k ≤ 1 2 3 3 Gradien t amplifiers, b fixed 2 2 – Cry ostat noise, small R 2 2 2 + k Cry ostat noise, large R 1 1 1 ∗ Larger in practice, b ecause of gradiometer im balance and field inhomogeneities. IV. SENSOR ARRA YS A. Combining data from m ultiple c hannels It is common to w ork with absolute v alues of the com- plex images to eliminate phase shifts. Images from mul- tiple c hannels can then be combined b y summing the squares and taking the square ro ot. This pro cedure, how- ev er, causes asymmetry in the noise distribution and loses information that can b e used for impro v ed com bination of the data. If the sensor arra y and the correlations of noise b et w een different sensors are known, the m ulti-channel data can b e com bined more effectively . In the following, we sho w that, where multiple sensors can form a softw are gradiometer, an array of N sensors can form an N th -order combination optimized to give the b est SNR for eac h v oxel. T o follo w the deriv ation in Ref. [39], consider a v oxel cen tered at ~ r , and N sensors indexed by j = 1 , 2 , ..., N . Based on Sec. I I A, each sensor has a unit magnetization image s j ( ~ r ) = β ∗ j ( ~ r ) V , where β j and V are the sensitiv- it y profile and vo xel v olume, respectively . Th e absolute v alue | s j | giv es the sensed signal amplitude caused by a unit magnetization in the vo xel, precessing p erpendicular to ~ B L . The complex phase represents the phase shift in the signal due to the geometry . T o study the p erformance of the array only , w e set V to unity . F or a vo xel cen tered at ~ r , we hav e a vector of recon- structed image v alues v = [ v 1 , v 2 , ..., v N ] > corresp onding to the N sensors. At this p oin t, the v alues v j ha ve not b een corrected according to the sensitivity . The linear com bination that determines the final v oxel v alue u can b e written in the form u = N X j =1 a ∗ j v j = a † v , (20) where † denotes the conjugate transpose. Requiring that the outcome is sensitivity-corrected sets a condition on the co efficient v ector a = [ a 1 , ..., a N ] > . In the absence of noise, a unit source magnetization gives v j = s j ( ~ r ) . The final v o xel v alue u should represent the source, whic h leads to the condition a † s = 1 . (21) Belo w, we sho w how a = [ a 1 , ..., a N ] > should b e chosen in order to maximize the SNR in the final image given the sensor array and noise prop erties. The single-sensor image v alues v i can b e written in the form v j = w j + ξ j where w j is the ‘pure’ signal and ξ j is the noise. The noise terms ξ j can b e modeled as random v ariables, which, for unbiased data, hav e zero exp ecta- tion: E( ξ j ) = 0 . If there is a bias, it can b e measured and subtracted from the signals b efore this step. The exp ectation of the final v alue of this vo xel is then E( u ) = E  a † ( w + ξ )  = a † w . (22) The noise in the vo xel is quantified by the v ariance of u . Eqs. (20) and (22) yield u = E( u ) + a † ξ , leading to V ar( u ) = E  | u − E( u ) | 2  = E h a † ξ ξ † a i = a † Σa , (23) where Σ = E( ξ ξ † ) identifies as the noise cov ariance ma- trix. F or simple cases, Σ is the same for all v oxels. How- ev er, it ma y v ary betw een vo xels if, for instance, the v ox- els are of different sizes. No w, the task is to minimize the noise a † Σa sub ject to the constraint in Eq. (21). The Lagrange multiplier metho d turns the problem into finding the minimum of L = a † Σa − λ (1 − a † s ) (24) with resp ect to a , while still requiring that Eq. (21) holds. F rom the constrain t it follows that a † s is real, so it ma y b e replaced by ( a † s + s † a ) / 2 in Eq. (24). By ‘completing the square’ in Eq. (24), one obtains L = ( a − ˜ a ) † Σ ( a − ˜ a ) − λ + constan t , (25) where ˜ a satisfies 2 Σ˜ a = − λ s . (26) Since Σ , b eing a co v ariance matrix, is p ositive (semi)definite, the minimum of L is found at a = ˜ a . F urther, Σ is alw ays inv ertible, as the contrary would imply that some non-trivial linear com bination of the signals would contain zero noise. Multiplying Eq. (26) b y s † Σ -1 from the left and using Eq. (21) leads to λ = − 2 / s † Σ -1 s . When this expression for λ is put bac k in to Eq. (26), the optimal choice for the coefficient vector a = ˜ a is obtained as a = Σ -1 s s † Σ -1 s . (27) 10 Similar to Eq. (7) of Ref. [56], Eqs. (23) and (27) reveal the final noise v ariance σ 2 fin for the given vo xel p osition, σ 2 fin = a † Σa = 1 s † Σ -1 s . (28) In the ab ov e deriv ation, we assumed little ab out how the individual single-sensor data w ere acquired. In fact, the only significant requiremen t was that the sensitivi- ties s i are well defined and accessible. As discussed pre- viously , the signal can b e mo deled to high accuracy at ULF (see Sec. I I A). B. Figures of merit and scaling for arrays Giv en the N th -order combination from Eqs. (20) and (27), the contribution of the sensor array to the v oxel- wise image SNR is giv en b y Eq. (28). W e define the arr ay-sensitivity-to-noise r atio aSNR as aSNR = √ s † Σ -1 s . (29) When each sensor in the arra y sees an equal flux noise lev el σ , the aSNR 1 / 2 tak es the form aSNR = √ s † X -1 s σ = arra y sensitivity noise level , (30) where X = Σ /σ 2 is the dimensionless noise c orr elation matrix. W e refer to the quan tity √ s † X -1 s as the ar- r ay sensitivity , which for weak correlation is given ap- pro ximately as || s || 2 . Scalin g law exp onen ts for the ar- ra y sensitivity are denoted by α a , and for the aSNR by δ = α a − α n . C. Correlation of noise b etw een sensors As already seen in Secs. IV A and IV B, the aSNR is affected by the correlation of random noise b etw een dif- feren t single-sensor channels. There are tw o main rea- sons for such correlations. First, a noise source that is not an intrinsic part of a sensor can directly couple to man y sensors. F or instance, thermal noise in conductors close to the sensors may result in such correlated noise (see Sec. II I B). Second, the pickups of the sensors them- selv es are coupled to each other through their m utual inductances. This cross-coupling increases noise correla- tion and ma y also affect the sensitivit y profiles via signal cross-talk. T o see the effect of noise correlation on the image SNR, consider a noise cov ariance matrix of the form Σ = σ 2 ( I + C ) , (31) where I is the identit y matrix and C contains the cor- relations b et ween channels (the off-diagonal elemen ts of X ). In words, each channel has a noise v ariance of σ 2 and c hannels p and q hav e correlation C pq = E( ξ p ξ ∗ q ) /σ 2 . Assume further that absolute v alues of the correlations C pq are substantially smaller than one. T o first order in C , the inv erse of Σ is obtained as Σ − 1 ≈ σ − 2 ( I − C ) . The SNR in the final image, accord- ing to Eq. (28), is then prop ortional to σ − 1 fin , with σ − 2 fin ≈ σ − 2  s † s − s † Cs  = σ − 2 k s k 2 2 − 2 σ − 2 X p 0 . This leads to the conclusion that the noise correlation tends to decrease the image SNR. While the assumptions made in the ab o ve discussion ma y not alwa ys b e exactly correct, the result is an in- dication that the correlation of noise b etw een adjacen t sensors is usually harmful—even if it is tak en in to ac- coun t in reconstruction. Moreov er, the actions tak en in order to reduce noise correlation are often such that the noise v ariances decrease as well. F or instance, eliminat- ing a noise source from the vicinity of the sensor array do es exactly that. Correlation can also b e reduced by minimizing the in ter-sensor cross-talk, for instance b y designing a sensor arra y with low mutual inductances betw een pic kup coils. If the mutual inductances are non-zero, the cross-talk can b e dramatically reduced by coupling the feedback of the SQUID flux-lo c ked lo op to the pickup circuit instead of more directly into the SQUID lo op [41]. This wa y , the sup ercurren t in the pickup coil stays close to zero at all times. In theory , the cross-talk of the flux signals can b e completely eliminated by this metho d. Correlated noise originating from sources far from the sub ject’s head and the sensor array can also b e atten- uated by signal processing metho ds prior to image re- construction. The signal sp ac e sep ar ation metho d (SSS) w as developed at Elekta Neuromag Oy [30] (now ME- GIN) for use with ‘whole-head’ MEG sensor arrays. The SSS metho d can distinguish betw een signals from inside the sensor helmet and those pro duced by distant sources. No w, the strong noise correlation is in fact exploited to significan tly impro ve the SNR. Similar methods may b e applicable to ULF MRI as w ell. T o help such metho ds, additional sensors can b e placed outside the helmet ar- rangemen t to provide an improv ed noise reference. F or sensor array comparisons, we assume that all mea- sures hav e b een taken to reduce correlated noise b efore image reconstruction. The details of the remaining noise correlation dep end on many , generally unknown aspects. Therefore, we set C = 0 in Eq. (31) for a sligh tly opti- mistic estimate, i.e. , sensor noises are uncorrelated, each ha ving v ariance σ 2 . 11 D. Filling the arra y In this section, w e use general scaling argumen ts to pro vide estimations of how the whole sensor array p er- forms as a function of the pickup coil size. Consider a surface, for instance, of the shap e of a helmet, and a vo xel at a distance l from the surface. The surface is filled with N pickup coils of radius R to measure the field perp en- dicular to the surface. W e assume the pickup coils are p ositioned either next to eac h other or in suc h a wa y that their areas ov erlap by a given fraction (see Fig. 1). The n umber of sensors that fit the surface is then prop ortional to R − 2 . T ake, at first, a v oxel far from the sensors; l  R . Now, the signal from the vo xel is spread ov er many sensors. F or Σ = σ 2 I , the aSNR is proportional to k s k 2 /σ . Assume that s j ∝ R α s and σ ∝ R α n , whic h leads to k s k 2 ∝ √ N R α s ∝ R α s − 1 , and finally , aSNR ∝ R δ , δ = α a − α n = α s − α n − 1 . (33) Here w e th us ha ve array sensitivity scaling according to α a = α s − 1 , as opp osed to α a = α s when N is fixed. Re- call from Sec. I I B that the flux sensitivities scale as R α s with α s = 2 for magnetometers and α s = 3 for first-order planar gradiometers, given that l  R . Assuming, for instance, optimally matched input circuits, the in trinsic flux noise of the sensor in b oth cases has a p o wer law b eha vior with exp onent α n = 1 / 2 (see Sec. II I B), which yields δ = 0 . 5 and δ = 1 . 5 . This is clearly in fa vor of using larger pickup coils. Especially for larger R , how- ev er, the cry ostat noise may b ecome dominan t, and one has α n ≈ 1 . Now, magnetometer arra ys hav e δ ≈ 0 , i.e. , the coils size do es not affect the SNR. Still, gradiometer arra ys p erform b etter with larger R ( α a ≈ 1 ). In the p erhaps unfortunate case that noise sources far from the sensors are dominan t, the noise behav es like the signal, that is, α s = α n and δ = − 1 . Unlike in the other cases, a higher SNR would be reached b y decreasing the pic kup coil size. Ho wev er, such noise conditions are not realistic in the lo w-correlation limit. Instead, one should aim to suppress the external noise b y improving the system design or by signal pro cessing. The breakdown of the assumption of l  R needs some atten tion. If the vo xel of interest is close to the sensor arra y , the image v alue is formed almost exclusively b y the closest pickup-loop. Now, for non-ov erlapping pickups, the results for single sensors ( α a = α s ) are applicable, and the optim um magnetometer size is R ≈ l . But then, if the vo xel is far from the array (deep in the head), and R is increased to the order of l , it is more difficult to dra w conclusions. W e therefore extend this discussion in Secs. V and VI by a computational study . V. METHODS FOR NUMERICAL STUD Y In order to b e able to c ompare the p erformance of dif- feren t sensor configurations, w e used 3-D computer mo d- els of sensor arrays and calculated their sensitivities to signals from different lo cations in the sample. The sensitivities of single pic kup coils were calculated using ~ B s from Eq. (8). Ev aluating the line integral re- quired the coil path ∂ S to be discretized. The n umber of discretization p oints could b e k ept small by analytically in tegrating Eq. (34) o ver n straigh t line segmen ts b e- t ween consecutive discretization points ~ r j and ~ r j +1 (the end p oin t ~ r n = ~ r 0 ): ~ B s ( ~ r ) = µ 4 π n − 1 X k =0 Z ~ r j +1 ~ r 0 = ~ r j d ~ r 0 × ( ~ r − ~ r 0 ) | ~ r − ~ r 0 | 3 . (34) As shown in App endix A, this integrates exactly to ~ B s ( ~ r ) = µ 4 π n − 1 X j =0 a j + a j +1 a j a j +1 ~ a j × ~ a j +1 a j a j +1 + ~ a j · ~ a j +1 , (35) where ~ a j = ~ r j − ~ r . Besides reducing computational com- plexit y and increasing accuracy , this result allo w ed exact computation for p olygonal coils. F or a precession field ~ B L = B L b e L , the single-sensor sensitivities w ere obtained from Eq. (10) and the array- sensitivit y and aSNR maps were computed according to Sec. IV B. The normalization of the v alues computed here is somewhat arbitrary; the real image SNR dep ends on a host of details that are not known at this p oin t (see Sec. I I A). How ever, the results can b e used for studying arra y sensitivity patterns and—with noise lev els scaled according to estimated coil inductances—for comparing differen t p ossible array setups. VI. RESUL TS Numerical calculations were p erformed for simple spherical sensor arra ys (Sec. VI A) as well as for real- istic configurations (Sec. VI B), e.g. , of the shap e of a helmet. The former were used for studying scaling b e- ha vior of arra y sensitivities with sensor size and n umber, extending the discussion in Sec. IV D. The latter w ere used for comparing array sensitivit y patterns of different p oten tial designs. A. Effects of size and num b er A sensor array mo del w as built b y filling the surface of a sphere of radius 10 cm (see Fig. 5) with N magnetome- ters or N / 2 planar units of tw o orthogonal planar first- order gradiometers. Combining one of the magnetome- ters with one of the gradiometer units would th us give a sensing unit similar to those of the Elekta/Neuromag MEG system, though circular (radius R ). All sensors w ere orien ted to measure the radial component of the field. A spherical surface of radius 6 cm was chosen to represen t the cerebral cortex. The cortex surface was 12 Figure 5. Geometry used in numerical analysis of the dep en- dence of array sensitivity as functions of sensor size R and n umber N at different p oin ts inside the imaging volume. Sen- sors are on a spherical surface of radius 10 cm. A shell with radius 6 cm is represen tativ e of points on the cerebral cortex. th us at distance 4 cm from the sensor shell. In addition, the center of the sphere was considered to represent deep parts of the brain. The data in Fig. 6 show the dep endence of the ar- ra y sensitivity on R . Note that the num b er of sensors is approximately prop ortional to R − 2 . The largest coil size R = 10 cm corresp onds to one magnetometer or gra- diometer unit on each of the six faces of a cub e. The solid lines corresp ond to the scaling of the sensitivit y as R α a , α a = α s − 1 . F or smaller R , the scaling laws from Sec. IV D hold in all cases, and particularly well for gra- diometers and deep sources. The scaling law fails most notably with the magnetometer array at the cortex. In- deed, the sensitivity starts to de cr e ase with R when R is v ery large, as was shown for a sp ecial case in Sec. I I B. The error bars in Fig. 6 corresp ond to the minimum and maximum v alue of the sensitivity at the cortex while the data symbols corresp ond to the a verage v alue. De- spite the strong orientational dep endence of single sen- sors (see Sec. II B), the arra y sensitivities are fairly uni- form at the cortex. Only at large R do the orien tational effects emerge. Figure 7 shows a different dataset on how the arra y sen- sitivit y changes with how densely the sensors are pack ed in to the array . In this case, a v arying num b er of mag- netometer coils or gradiometer units with fixed radius R = 1 . 44 cm w as distributed on the spherical shell. The aSNR of v oxels at the cen ter scales as √ N to an excel- len t accuracy . While the av erage sensitivity at p oin ts on the cortex also ob eys √ N scaling remark ably w ell, the uniformit y drops dramatically when N is low ered below roughly 30 sensors. Closer to the sensors, e.g. on the scalp, this effect is even more pronounced. 10 − 2 10 − 1 Radius R of pic kup coils [m] 10 0 10 1 10 2 Arra y sensitivit y or noise lev el α n = 1 α n = 1 2 α a = 2 α a = 1 α a = 3 Cortex, magn. Cen ter, magn. Cortex, grad. Cen ter, grad. Figure 6. Scaling of arra y sensiti vit y at the center and on the cortex as depicted in Fig. 5: sphere filled with magnetometer lo ops and with planar units of t wo orthogonal gradiometers arranged side by side. Error bars corresp ond to the minimum and maxim um v alues. Noise scaling with size is included in the figure, illustrating a p otential cross-o ver from sensor noise with α n = 1 / 2 to cryostat noise or suboptimal input circuit matc hing with α n = 1 . With fixed N , array sensitivit y scaling is steep er and giv en by α a = 2 , 3 for planar magnetometers and gradiometers. 10 1 10 2 No. magnetometers N / 1 10 0 10 1 Arra y sensitivit y Cortex Cen ter 10 1 10 2 No. gradiometers N / 2 10 0 10 1 ∝ √ N Figure 7. Scaling of array sensitivity as √ N at the center and on the cortex as depicted in Fig. 5, when the pickup coil ra- dius is fixed at R = 1 . 44 cm: (left) N magnetometers, (right) N/ 2 planar units of tw o orthogonal gradiometers. Error bars corresp ond to the minimum and maximum v alues. B. Realistic sensor configurations Figure 8 presents several p ossible sensor configurations and pro vide maps of log 10 ( aSNR ) for their comparison for their comparison. The data sho wn are sagittal slices of the 3-D maps, i.e. , on the symmetry plane of the sensor 13 arra y . Other slices, ho wev er, displa yed similar perfor- mance at the cortex. Also changing the direction of the precession field ~ B L had only a minor effect on the SNR in the region of in terest. In all cases shown here, ~ B L w as parallel to the y axis, which is p erp endicular to the visualization plane. Note that this contrasts MRI con- v ention, where the ~ B L direction is considered fixed and alw ays along the z axis. In most cases, the sensors are arranged on a helmet sur- face at 102 p ositions as in the Elekta/Neuromag system. Again, magnetometers and planar double-gradiometer units are considered separately (here, R = 1 . 25 cm, re- sem bling conv entional MEG sensors). The same flux noise level was assumed for magnetometers and planar gradiometers of the same size. In addition, we consider arra ys with axial gradiometers as well as radially oriented planar gradiometers, b oth cases having k = 1 , b = 4 cm and R = 1 . 25 cm. Configurations with 102 ov erlapping units with R = 2 . 5 cm are also considered, as well as the existing Los Alamos 7-channel coil geometry [42] and the single large second-order gradiometer at UC Berkeley [2] (see figure caption). F or long-baseline gradiometers with k = 1 , L p w as estimated to b e twice that of a single loop, and six times for k = 2 . With planar sensor units of R = 1 . 25 cm [Fig. 8(a–b)], the aSNR for 102 magnetometers is three times that of 204 gradiometers at the cerebral cortex. At the center of the head, the difference is almost a whole order of magnitude in fa vor of the magnetometers. Therefore, the small gradiometers bring little improv ement to the image SNR if the magnetometers are in use. How ever, as shown previously , esp ecially gradiometer p erformance impro ves steeply with coil size. Allo wing the coils to o verlap with R = 2 . 5 cm [Fig. 8(g–h)] leads to a v astly impro ved aSNR, esp ecially with gradiometers, but also with magnetometers. Gradiometers with long baselines provide somewhat magnetometer-lik e sensitivity patterns while rejecting ex- ternal noise. How ev er, their aSNR p erformance is infe- rior to magnetetometers b ecause of their larger induc- tance, yielding higher flux noise when the sensor noise dominates; see Sec. I I I A. Helmet arrays of magnetome- ters can provide a similar aSNR in the deep est parts of the brain as the Berkeley gradiometer pro vides at a small area on the scalp. VI I. CONCLUSIONS AND OUTLOOK Extending Ref. [39], we analyzed a v ariety of factors that affect the noise and sensitivit y of a SQUID-based sensor array for ULF MRI of the brain. Many of the principles, how ever, apply to non-SQUID arra ys as well. W e also deriv ed numerical means for studying and com- paring the SNR p erformances of an y given sensor array designs. Signal- and noise-scaling arguments and calculations sho wed that filling a sensor array with a huge n um b er of tin y sensors is usually not adv antageous. Larger pickup coil sizes give a b etter image SNR at the center of the head and, up to some p oin t, also at closer sources such as the cerebral cortex. This is true even if the n umber of sensors needs to b e decreased due to the limited area a v ailable for the arra y . How ever, the av erage v o xel SNR is prop ortional to the square root of the num b er of sensors. Sev eral p ossible arra y designs w ere compared, includ- ing existing arrays designed for MEG and ULF MRI. The results are mostly in fav or of magnetometers and large first-order gradiometers. While t ypically having inferior SNR, gradiometers do hav e the adv antage of rejecting ex- ternal fields, reducing also transien t issues due to pulsed fields [52]. An esp ecially dramatic difference was found when comparing a magnetometer-filled helmet with a sin- gle larger gradiometer. In general, using an arra y of sensors relaxes the dy- namic range requirements for sensor readout. Splitting a large lo op into smaller ones further allo ws interference rejection based on correlation, while also increasing the SNR close to the center of the lo op. An array of many sensors also solv es the single-sensor problem of ‘blind an- gles’. Our initial analysis of overlapping magnetometer and gradiometer coils gav e promising results. Implementing suc h arra ys, how ev er, p oses challenges. Practical con- siderations include ho w to fabricate suc h an array and what materials to use. F or instance, wire-wound Type-I sup erconducting pic kup coils hav e sho wn some fav orable prop erties [57, 58] in pulsed systems, and exploiting the dynamics of sup erconductor-p enetrating flux [52, 53, 59] has b een promising. How ever, existing techniques are not suitable for helmet configurations with o v erlapping coils. In addition, careful design w ork should b e conducted to minimize m utual inductances and other coupling issues. F urther significant improv emen ts could b e achiev ed by placing the sensors closer to the scalp, but that w ould require dramatic adv ancemen ts in cryostat tec hnology , and was not studied here. Here, we only considered the contribution of the sensor arra y to the imaging performance. Other things to con- sider are the p olarizing technique as well as the ability of the instrumentation to apply more sophisticated se- quences and reconstruction techniques, while preserving lo w system noise. A class of techniques enabled by m ulti- c hannel magnetometers is accelerated parallel MRI [32]. Ho wev er, the so-called geometry factor should b e taken in to account [60] if large parallel acceleration factors are pursued. [1] A. Abragam, The principles of nucle ar magnetism , 32 (Oxford Universit y Press, New Y ork, USA, 1961). [2] J. Clarke, M. Hatridge, and M. Mößle, Annual Reviews 14 Figure 8. Base-10 logarithms of aSNR for differen t sensor-array geometries. T o allo w comparison of different arra ys, we assumed SQUID noise scaling according to optimally matched input circuits. (a) Magnetometers: R = 1 . 25 cm, (b) double-gradiometer units: R = 1 . 25 cm, (c) axial gradiometers: b = 4 cm, R = 1 . 25 cm, (d) 7 Los Alamos second-order axial gradiometers: b = 6 cm, R = 1 . 85 cm, (e) Berkeley single second-order axial gradiometer: b = 7 . 5 cm, R = 3 . 15 cm, (f ) radially oriented planar gradiometers [Fig. 2(f )]: b = 4 cm, R = 1 . 25 cm, (g) ov erlapping double-gradiometer units: R = 2 . 5 cm, (h) ov erlapping magnetometers: R = 2 . 5 cm. The data rate of the acquisition is prop ortional to the square of the of aSNR. 15 in Biomedical Engineering 9 , 389 (2007). [3] A. Gramfort, M. Luessi, E. Larson, D. A. Engemann, D. Strohmeier, C. Bro db ec k, L. Parkk onen, and M. S. Hämäläinen, NeuroImage 86 , 446 (2014). [4] Y.-C. Hsu, K. C. J. Zev enhov en, Y.-H. Ch u, J. Dab ek, R. J. Ilmoniemi, and F.-H. Lin, Magnetic Resonance in Medicine 75 , 2255 (2016). [5] P . T. V esanen, K. C. J. Zev enhov en, J. O. Nieminen, J. Dabek, L. T. P arkkonen, and R. J. Ilmoniemi, J. Magn. Reson. 235 , 50 (2013). [6] S.-J. Lee, K. Kim, C. S. Kang, S. min Hwang, and Y.-H. Lee, Sup erconductor Science and T echnology 23 , 115008 (2010). [7] K. Buck enmaier, M. Rudolph, P . F ehling, T. Steffen, C. Bac k, R. Bernard, R. Pohmann, J. Bernarding, R. Kleiner, D. Koelle, et al. , Review of Scientific Instru- men ts 89 , 125103 (2018). [8] P . T. V esanen, J. O. Nieminen, K. C. Zev enho ven, Y.-C. Hsu, and R. J. Ilmoniemi, Magnetic Resonance Imaging 32 , 766 (2014). [9] J. O. Nieminen, K. C. J. Zevenho ven, P . T. V esanen, Y. C. Hsu, and R. J. Ilmoniemi, Magn. Reson. Imaging 32 , 54 (2014). [10] P . Hömmen, J.-H. Storm, N. Höfner, and R. Körber, Magnetic Resonance Imaging 60 , 137 (2019). [11] S. K. Lee, M. Mößle, W. Myers, N. Kelso, A. H. T rab esinger, A. Pines, and J. Clarke, Magnetic Reso- nance in Medicine 53 , 9 (2005). [12] S. Hart wig, J. V oigt, H.-J. Scheer, H.-H. Albrech t, M. Burghoff, and L. T rahms, The Journal of Chemical Ph ysics 135 , 054201 (2011). [13] P . T. V esanen, K. C. J. Zev enhov en, J. O. Nieminen, J. Dab ek, L. T. Parkk onen, and R. J. Ilmoniemi, Journal of Magnetic Resonance 235 , 50 (2013). [14] R. H. Kraus Jr., P . V olego v, A. Matlacho v, and M. Espy , NeuroImage 39 , 310 (2008). [15] R. K örb er, J. O. Nieminen, N. Höfner, V. Jazbinšek, H.- J. Sc heer, K. Kim, and M. Burghoff, Journal of Magnetic Resonance 237 , 182 (2013). [16] Y. Xue, J.-H. Gao, and J. Xiong, NeuroImage 31 , 550 (2006). [17] K. Kim, S.-J. Lee, C. S. Kang, S.-M. Hwang, Y.-H. Lee, and K.-K. Y u, NeuroImage 91 , 63 (2014). [18] K. Kim, AIP Adv ances 2 , 022156 (2012). [19] O. V. Lounasmaa and H. Seppä, Journal of Low T emp er- ature Physics 135 , 295 (2004). [20] J. V rba and S. E. Robinson, Sup erconductor Science T echnology 15 , R51 (2002). [21] V. Pizzella, S. Della P enna, C. Del Gratta, and G. L. Ro- mani, Sup erconductor Science and T echnology 14 , R79 (2001). [22] M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Kn uutila, and O. V. Lounasmaa, Reviews of Modern Ph ysics 65 , 413 (1993). [23] C. Del Gratta, V. Pizzella, F. T ecchio, and G. L. Ro- mani, Rep orts on Progress in Physics 64 , 1759 (2001). [24] P . T. V esanen, J. O. Nieminen, K. C. J. Zevenho ven, J. Dab ek, L. T. Parkk onen, A. V. Zhdanov, J. Luoma- haara, J. Hassel, J. Pen ttilä, J. Simola, et al. , Magnetic Resonance in Medicine 69 , 1795 (2013). [25] P . E. Magnelind, J. J. Gomez, A. N. Matlasho v, T. Owens, J. H. Sandin, P . L. V olegov, and M. A. Esp y , IEEE T ransactions on Applied Sup erconductivity 21 , 456 (2011). [26] J. Luomahaara, M. Kiviranta, L. Grön b erg, K. C. J. Zev enhov en, and P . Laine, IEEE T ransactions on Ap- plied Sup erconductivit y 28 , 1600204 (2018). [27] R. Körber, J.-H. Storm, H. Seton, J. P . Mäkelä, R. Pae- tau, L. P arkkonen, C. Pfeiffer, B. Riaz, J. F. Schneider- man, H. Dong, et al. , Sup erconductor Science and T ech- nology 29 , 113001 (2016). [28] A. J. Mäkinen, K. C. Zev enhov en, and R. J. Ilmoniemi, IEEE T ransactions on Medical Imaging 38 , 1317 (2019). [29] M. A. Uusitalo and R. J. Ilmoniemi, Medical and Biolog- ical Engineering and Computing 35 , 135 (1997). [30] S. T aulu and M. Ka jola, Journal of Applied Physics 97 , 124905 (2005). [31] K. P . Pruessmann, M. W eiger, M. B. Scheidegger, and P . Bo esiger, Magnetic Resonance in Medicine 42 , 952 (1999). [32] D. J. Larkman and R. G. Nunes, Physics in Medicine and Biology 52 , R15 (2007). [33] M. A. Ohliger and D. K. Sodickson, NMR in Biomedicine 19 , 300 (2006). [34] K. C. J. Zevenho ven and S. Alanko, Journal of Ph ysics: Conference Series 507 , 042050 (2014). [35] A. I. Ahonen, M. S. Hämäläinen, R. J. Ilmoniemi, M. J. Ka jola, J. E. T. Knuutila, J. T. Simola, and V. A. Vilk- man, IEEE T ransactions on Biomedical Engineering 40 , 859 (1993). [36] J. Nurminen, The magnetostatic multip ole exp ansion in biomagnetism: applic ations and implic ations , Ph.D. the- sis, Aalto Universit y , Finland (2014). [37] E. Burmistrov, A. Matlasho v, H. Sandin, L. Sc hultz, P . V olegov, and M. Espy , IEEE T ransactions on Applied Sup erconductivit y 23 , 1601304 (2013). [38] W. Myers, D. Slich ter, M. Hatridge, S. Busc h, M. Mößle, R. McDermott, A. T rab esinger, and J. Clark e, Journal of Magnetic Resonance 186 , 182 (2007). [39] K. Zevenho ven and R. J. Ilmoniemi, Pro ceedings of the In ternational So ciet y of Magnetic Resonance in Medicine 19 , 1803 (2011). [40] Y.-C. Hsu, P . T. V esanen, J. O. Nieminen, K. C. J. Zev enhov en, J. Dab ek, L. Parkk onen, I.-L. Chern, R. J. Ilmoniemi, and F.-H. Lin, Magnetic Resonance in Medicine 71 , 955 (2014). [41] J. Clarke and A. I. Braginski (Eds.), The SQUID Hand- b o ok (Wiley-VCH V erlag Gm bH & Co. KGaA, W einheim, German y , 2004). [42] V. S. Zotev, A. N. Matlac ho v, P . L. V olegov, H. J. Sandin, M. A. Espy , J. C. Mosher, A. V. Urbaitis, S. G. Newman, and R. H. Kraus, Jr., IEEE T ransactions on Applied Su- p erconductivit y 17 , 839 (2007). [43] J. Mates, K. Irwin, L. V ale, G. Hilton, and H. Cho, Journal of Low T emp erature Physics 176 , 483 (2014). [44] F. W. Grov er, Inductanc e Calculations: W orking for- mulas and tables (Dov er Publications, Inc., New Y ork, 1973). [45] H. Nyquist, Physical Review 32 , 110 (1928). [46] J. B. Johnson, Physical Review 32 , 97 (1928). [47] J. Storm, P . Hömmen, N. Höfner, and R. Körber, Mea- suremen t Science and T echnology 30 , 125103 (2019). [48] H. C. Seton, J. M. S. Hutchison, and D. M. Bussell, Cry ogenics 45 , 348 (2005). [49] A. T ervo, Noise optimization of multi-layer insulation in liquid-helium cryostat for br ain imaging , Master’s thesis, Aalto Universit y , Esp o o, Finland (2016). [50] T. V arpula and T. P outanen, Journal of Applied Physics 16 55 , 4015 (1984). [51] J. Nenonen, J. Mon tonen, and T. Katila, Review of Sci- en tific Instruments 67 , 2397 (1996). [52] K. C. J. Zevenho ven, Solving tr ansient pr oblems in ultr a- low-field MRI , Master’s thesis, Universit y of California, Berk eley and Aalto Univ ersity , Finland (2011). [53] K. C. J. Zevenho ven et al. , “ T ow ards high-quality ultra- lo w-field MRI with a sup erconducting polarizing coil,” P ap er 3A-EL-O4, presented at the 11th Eur. Conf. Appl. Sup ercond., Geno v a, Italy . (2013). [54] K. C. J. Zev enhov en, S. Busch, M. Hatridge, F. Öisjöen, R. J. Ilmoniemi, and J. Clark e, Journal of Applied Ph ysics 115 , 103902 (2014). [55] K. C. J. Zev enhov en, H. Dong, R. J. Ilmoniemi, and J. Clarke, Applied Ph ysics Letters 106 , 034101 (2015). [56] J. Cap on and N. Go odman, Pro ceedings of the IEEE 58 , 1785 (1970). [57] J. Luomahaara, P . V esanen, J. P enttilä, J. Nieminen, J. Dab ek, J. Simola, M. Kiviranta, L. Grönberg, C. J. Zev enhov en, R. Ilmoniemi, et al. , Superconductor Sci- ence and T echnology 24 , 075020 (2011). [58] S.-m. Hwang, K. Kim, K. Kyu Y u, S.-J. Lee, J. Hyun Shim, R. Körber, and M. Burghoff, Applied Ph ysics Letters 104 , 062602 (2014). [59] E. Al-Dabbagh, J.-H. Storm, and R. Körber, IEEE T ransactions on Applied Sup erconductivit y 28 , 1 (2018). [60] F.-H. Lin, P . T. V esanen, J. O. Nieminen, Y.-C. Hsu, K. C. J. Zevenho v en, J. Dabek, L. T. P arkkonen, A. Zh- dano v, and R. J. Ilmoniemi, Magnetic Resonance in Medicine 70 , 595 (2013). App endix A: Exact Biot–Sa v art in tegral o v er p olygonal path Here, w e derive an exact expression for calculating the Biot–Sa v art in tegral ov er a p olyline, i.e. a path consisting of connected line segments; see Eq. (34). Consider a line segmen t from ~ r p to ~ r q . Using the notations ~ a = ~ r − ~ r 0 and ~ a j = ~ r − ~ r j ( j = p, q ), the integrals in Eq. (34) can b e written as I = Z ~ r q ~ r 0 = ~ r p d ~ r 0 × ~ a a 3 = Z ~ a q ~ a = ~ a p ~ a × d ~ a a 3 = Z 1 0 [ ~ a p + ( ~ a q − ~ a p ) t ] × ( ~ a q − ~ a p ) | ~ a p + ( ~ a q − ~ a p ) t | 3 dt , (A1) where the line segmen t has b een parametrized as ~ a = ~ a p + ( ~ a q − ~ a p ) t , t ∈ [0 , 1] . F or an arbitrary vector ~ V , one has ~ V × ~ V = 0 . Apply- ing this twice to Eq. (A1) (with ~ V = ~ a q − ~ a p and ~ V = ~ a p ) yields I = Z 1 0 ~ a p × ~ a q | ~ a p + ( ~ a q − ~ a p ) t | 3 dt . (A2) The divisor can b e expanded as { [ ~ a p + ( ~ a q − ~ a p ) t ] 2 } 3 2 = ( C t 2 + Dt + E ) 3 2 , (A3) with      C = ( ~ a q − ~ a p ) 2 , D = 2 ~ a p · ( ~ a q − ~ a p ) , E = ~ a 2 p . (A4) The relev ant integral is given b y ˜ I = Z 1 0 ( C t 2 + Dt + E ) − 3 2 dt =  2(2 C t + D ) (4 C E − D 2 ) √ C t 2 + Dt + E  1 t =0 , (A5) as can b e verified b y differen tiation. Straigh tforward al- gebraic manipulation leads to simplified expressions:          4 C E − D 2 = 4 a 2 p a 2 q − 4( ~ a p · ~ a q ) 2 , 2 C + D = 2 ~ a q · ( ~ a q − ~ a p ) , √ C + D + E = a q , √ E = a p . (A6) No w, the integral in Eq. (A5) b ecomes ˜ I = 1 a 2 p a 2 q − ( ~ a p · ~ a q ) 2  ~ a q · ( ~ a q − ~ a p ) a q − ~ a p · ( ~ a q − ~ a p ) a p  , whic h simplifies as follows: ˜ I = 1 a p a q ( a p ~ a q − a q ~ a p ) · ( ~ a q − ~ a p ) a 2 p a 2 q − ( ~ a p · ~ a q ) 2 = 1 a p a q ( a p + a q )( a p a q − ~ a p · ~ a q ) ( a p a q + ~ a p · ~ a q )( a p a q − ~ a p · ~ a q ) = 1 a p a q a p + a q a p a q + ~ a p · ~ a q . (A7) Using the final expression in Eq. (A7), the original in- tegral I = ˜ I ( ~ a p × ~ a q ) can b e written as I = a p + a q a p a q ~ a p × ~ a q a p a q + ~ a p · ~ a q , (A8) pro ving the identit y of Eq. (35). In addition to calcu- lating exact Biot–Sav art integrals for p olylines, Eq. (35) can also b e used for efficient numerical integration o ver arbitrary discretized paths. A Python pack age emfields optimized for efficient computation will b e released on the Python Pac k age Index (PyPI).