In MRI, the Larmor precession frequency ω(x, y, z) = γB(x, y, z) of the proton spins in the position-dependent magnetic field B(x, y, z) frequency-and phase-encodes the proton density distribution into a magnetic signal that is subsequently decoded to form an image 4 (γ is the magnetogyric ratio). In clinical MRI machines 4 the strength of the applied homogeneous static magnetic field B 0 = B 0 ẑ is typically 1.5 T. There has been recent interest, however, in systems operating in magnetic fields of the order of 10 -4 T (for example [13][14][15][16][17][18] ), where T 1 -weighted contrast is significantly enhanced 16 (T 1 is the longitudinal relaxation time). The loss of polarization is compensated-at least in part-by 3 prepolarizing 19 the spins at a much higher field, or by hyperpolarization techniques using lasers 20 , dynamic nuclear polarization 21,22 or parahydrogen-induced polarization 23 . The loss of signal amplitude inherent in Faraday-Law detection is mitigated by detecting the nuclear magnetization with either a Superconducting QUantum Interference Device (SQUID) 9 or an atomic magnetometer 24 , both of which respond to the magnetic flux itself, rather than its time rate of change. Regardless of the magnitude of B 0 , all currently used imaging processes involve the superposition of magnetic field gradients on a static field to impose spatial variations of the total field across the subject or sample. In the zero static field regime reported here, conventional MRI gradients are unable to encode the spins along a given direction and Fourier encoding breaks down.In conventional MRI techniques, the applied magnetic field gradients are assumed to be linear and unidirectional so that the field due to gradients is given by B(x, y, z) = (G x x + G y y + G z z) ẑ , where G x = ∂B z /∂x, G y = ∂B z /∂y, and G z = ∂B z /∂z are constants 4 . As an example, B(x, y, z) = G z z ẑ is shown in Fig. 1a. In reality, however, such idealized gradients are forbidden by the Maxwell equations divB = curlB = 0 for any magnetic field B in free space. In fact, any gradient must be accompanied by concomitant gradients in at least one additional direction, as illustrated in Fig The gradient field B(y, z) = (∂B y /∂y)y ŷ + (∂B z /∂z)z ẑ is turned on, and subsequently turned off nonadiabatically at time τ (point B). During this time interval, the spin precesses about B(y′, z′). The time τ is chosen to satisfy the requirement τ << 1/γG z L.Consequently, the precession during the interval τ is small, and we can treat it as the sum of precessions around ẑ and ŷ : δ z = γ(∂B z /∂z)z'τ around ẑ (Fig. 2b) and δ y = γ(∂B y /∂y)y'τ around ŷ (Fig. 2c). After the gradient pulse, a π pulse of uniform field Β π is applied along the z-axis with amplitude and duration adjusted to produce a precession angle of π around ẑ . This pulse flips the spin to the point C in Figs. 2b and 2c.
5Subsequently, a second gradient pulse brings the spin to D, and a second π pulse to E. This sequence of pulses produces a net precession of the spin about B z , ...