key: cord-0221224-hglju9y6
authors: Lehmkuhl, Soren; Fleischer, Simon; Lohmann, Lars; Rosen, Matthew S.; Chekmenev, Eduard Y.; Adams, Alina; Theis, Thomas; Technology, Stephan Appelt Institute of Microstructure; Technology, Karlsruhe Institute of; Germany,; Chemistry, Department of; University, North Carolina State; Raleigh,; NC,; Technical, Institute of; Chemistry, Macromolecular; University, RWTH Aachen; Hospital, Massachusetts General; Imaging, A. A. Martinos Center for Biomedical; Boston,; MA,; USA,; Physics, Department of; University, Harvard; Cambridge,; Biosciences, Integrative; Institute, Karmanos Cancer; Detroit, Wayne State University; MI,; Sciences, Russian Academy of; Moscow,; Russia,; Raleigh, North Carolina State University; Engineering, Joint Department of Biomedical; Hill, University of North Carolina at Chapel; Engineering, Central Institute for; Electronics,; Systems, Analytics Electronic; GmbH, Forschungszentrum Julich; Julich,
title: RASER MRI: Magnetic Resonance Images formed Spontaneously exploiting Cooperative Nonlinear Interaction
date: 2022-03-01
journal: nan
DOI: nan
sha: 5a311232fcc5b62639fc717da35d9e83aa290570
doc_id: 221224
cord_uid: hglju9y6

The spatial resolution of magnetic resonance imaging (MRI) is fundamentally limited by the width of Lorentzian point spread functions (PSF) associated with the exponential decay rate of transverse magnetization (1/T2*). Here we show a different contrast mechanism in MRI by establishing RASER (Radio-frequency Amplification by Stimulated Emission of Radiation) in imaged media. RASER imaging bursts emerge out of noise and without applying (Radio Frequency) RF pulses when placing spins with sufficient population inversion in a weak magnetic field gradient. A small difference in initial population inversion density creates a stronger image contrast than conventional MRI. This contrast is based on the cooperative nonlinear interaction between all slices. On the other hand, the cooperative nonlinear interaction gives rise to imaging artifacts, such as amplitude distortions and side lobes outside of the imaging domain. Both the contrast and the artifacts are demonstrated experimentally and predicted by simulations based on a proposed theory. This theory of RASER MRI is strongly connected to many other distinct fields related to synergetics and non-linear dynamics.

includes inverted states (positive d 0 ), and it has been studied extensively in NMR spectroscopy (24, 27, (29) (30) (31) .

In order to understand how the RASER can be utilized for MR imaging, we introduce an analysis of the RASER action in the presence of a magnetic field gradient G z . The gradient creates a frequency range Δ = γ H ·G z · L that spans the image domain of the object of length L (SI, section 1). The initial nuclear spin population inversion is spread over the imaging domain Δ and is

, where (ν) is the population inversion density and ν 0 is the offresonance frequency in the center of the imaging domain. The integrand (ν)dν can be described as the number of negatively polarized spins in the frequency interval [ν,ν+dν] . Given a profile ( ), a total RASER MRI signal emerges spontaneously out of the spin noise. To generate a system where a numerical evaluation is feasible, we divide the image domain  into N = / individual slices. To avoid numerical artifacts, the distance δν between consecutive slices has to be chosen small enough. Specifically δν < w has to be fulfilled, where w = 1/(πT 2 * ) is the natural linewidth. Furthermore, to estimate whether a given d 0 is RASER active in a given gradient G z , we also introduce the threshold population density ρ d th = d th /w as used below. To calculate the dynamics of the nonlinearly coupled slices, each slice μ=1…N is characterized (4) .

The coupling constant is given as = 0 ℏ H 2 /(4V s ). The model for RASER MRI represented by Eqs.(1-4) is formulated in the rotating frame (for a complete derivation see SI, section 1), and is a modification of the existing multi-mode RASER theory (12, 16) . The modifications comprise the initial boundary conditions for μ (0) in Eq.(4), the absence of pumping in Eq.(1) and the definition of the slice frequencies in Eq. (3) . We assume a random fluctuation (nuclear spin noise excitation) for the initial small amplitudes A μ (0) and phases ϕ μ (0), which initiate the self-induced RASER burst in the absence of any RF excitation. Numerical simulations of Eqs.(1-4) reveal three important invariance principles for RASER MRI: Provided that δν < w and T 1 >> T 2 * , the amplitude and contrast properties of the RASER images are 4 independent of (I) the value of the slicing δν, (II) the longitudinal relaxation time T 1 and (III) the values of the initial conditions A  (0) and   (0).

Certain processes can be identified by examining the dynamics described by Eqs (1-3): The population inversion of a given mode μ in Eq.

(1) decays with the rate 1/T 1 and is decreased by the rate given by the sum over all quadratic terms −4 cos( − ). In turn, the amplitude of A μ in Eq.(2) decays with the rate 1/T 2 * and increases for τ = μ with the rate β·d μ . The last term on the right side of Eq.(2), μ ∑ τ cos( τ − μ ) τ=1 for τ ≠ μ involves a sum over all other amplitudes τ cos( τ − μ ). This sum can be a growth or decay rate for A μ , depending on the specific values of all other phase differences τ − μ . The collective action of all modes strongly influences the amplitude and sign of the rate dA μ /dt, which defines the amplitudes A μ of the final image.

The spatial encoding of each slice μ = 1..N is reflected by the first term in Eq. (3) , where each slice is oscillating at the angular frequency ω μ = 2π(ν 0 -0.5(Δ-δν(2μ-1))). Apart from this linear evolution of ϕ μ with time t, there is a nonlinear collective term ( μ / μ ) ∑ τ sin ( τ − μ ) =1 which is responsible for synchronism. In fact Eq. (3) is analogue to Kuramoto's model of synchronized oscillators (32-34). The dynamics of RASER MRI given by Eqs.

(1-4) can be described by a collection of synchronized oscillators or slices with distinct angular frequencies ω μ , where the amplitude A μ of each oscillator depends on the self-organization controlled by the collective interaction with all other slices. Therefore, the derivative of the amplitude of each slice depends on the mean-field amplitude produced by all other slices. Finally, the total RASER signal is obtained by the sum of all transverse spin components

. Here, we focus on the difference between the concept of single PSF´s to analyze conventional MR image formation. and the collective meanfield approach, which is the basis of RASER MRI. Numerical solutions of Eqs.(1-4) are evaluated (see Fig.1 ) in order to highlight the difference of the spin dynamics for a single RASER slice and the collective behavior of coupled slices.

The simplest case is shown in Fig. 1(A) for N = 1 and T 1 = , where the numerically evaluated form matches the exact solution introduced by Mao et al. (29, 35, 36) also discussed by others (37, 38) . The corresponding phased and absolute spectrum of α = α 1 are displayed on the bottom right of Fig.1(A) . For this case T 1 = , the PSF is a hyperbolic secant with width w sech (SI, section 2, Eq. S19). Close to the threshold, such a PSF is narrower than the Lorentzian NMR linewidth w = 1/(πT 2 * ), because the RASER signal involves dedamping.

No exact solution exists for a finite T 1 , but the MR signal represents an asymmetrically shaped PSF ( Fig. 1(B) , see SI, section 3). The linewidth w as in the spectrum is slightly broader compared to the symmetric case ( Fig.1(A) ), but still smaller than w. Here, for the first time, we include both the effects of finite T 1 and the nonlinear interactions between N slices formed in the presence of a gradient. In contrast to standard MRI, the image contrast and the spatial resolution cannot be explained by independent individual PSFs. Each slice is sensitive to the collective action of all other slices, which makes RASER imaging highly sensitive to local variations in d μ [SI, section 4(a)] providing new frontiers in MRI spin physics. 5 

In the simulation in Fig. 1(C) , a rectangular polarization profile (inset upper right) is assumed to generate a RASER signal in the presence of a field gradient. The time evolution of three of the N = 30 slices is depicted on the left. The shape of signal of these slices differs significantly from the uncoupled PSFs in (A) and (B). A corresponding 1D RASER image (projection) is obtained as the Fourier transform from Sig(t) = Re(∑α μ ). The amplitude in the center of the RASER image is larger and decaying side lobes arise outside of the image boundaries at x = 4 mm (bottom right). These artifacts are expected from the theory described in Eqs.(1-4) and evaluated in detail by numerical simulations in the SI, section 4. In Fig. 1(D) , we simulate a RASER image using a spin density profile ρ d (ν) to match the experimental setup described below in Fig. 2(A,B) . This non-uniform spin density profile ρ d (ν) entails two equal compartments separated by a gap. The evolution of five representative RASER slices out of N = 50 coupled slices is shown ( Fig. 1(D) , left). The image after Fourier transformation (bottom right) reflects roughly the shape ρ d (ν) except for the deformed amplitudes of the flat tops and the side lobes, which occur outside the imaging boundaries. The nonlinear interaction between all slices is mediated by the virtual photons (wavy arrows, wavelength ≫ sample dimension) in the resonator (red) (39). After the RASER burst, the Zeeman energy of the spins is fully transferred to the current of the coil (40). (A) For N = 1 and T 1 = , the signal α 1 (t) = α(t) is plotted in the (t, Re(α), Im(α)) space (left). The projection Re(α) for d 0 = 4.210 15 and the corresponding Fourier transformed spectra are shown on the right. (B) For N = 1, T 1 = 5 s and d 0 = 4.210 15 , the signal burst α(t) is asymmetric with respect to time. (C) Sketch of three representative signals α μ , μ = 1, 15, 30 out of N = 30 interacting slices (T 1 = 5 s, Δ = 6 Hz, rectangular profile with ρ d (ν) = 7.510 15 /Hz). (D) Five representative signals α μ out of N = 50 coupled slices (T 1 = 5 s, Δ = 10 Hz, non-uniform density ρ d (ν)). Threshold population density ρ d th = d th /w = 6.6 10 15 /Hz is indicated as dotted line in the insets in (C,D).

To experimentally examine the RASER MRI theory, a simple phantom was prepared consisting of a cylindrical sample chamber divided into two measurement chambers by a glass slide ( Fig.  2(A,B) ). The two chambers are individually supplied with p-H 2 to generate highly negative polarized proton spins (i.e. d 0 >> d th ). The chemical system chosen is pyrazine in a liquid methanol-d 4 solution with a dissolved iridium-based SABRE catalyst for nuclear spin polarization (18, 41). RASER MR images were acquired in the presence of weak G x and G z magnetic field gradients on the order of a few mG/cm. Conventional MR images (SEI) were obtained with a 90°-180° Spin Echo sequence ( Fig. 2(C) ) as a reference. Prior to the acquisition of the reference Spin Echo Image (SEI), a crusher field gradient was applied to the hyperpolarized sample, to suppress spontaneous RASER build-up. 1D images were acquired using the G z gradient in order to visualize the two chambers separated by the dividing glass slide. 2D images were recorded through stepwise switching of the G x and G z gradients to rotate through a circle with constant absolute gradient (G =(G z 2 +G x 2 ) 1/2 . The 2D image was then obtained via projection reconstruction, which is also common in Computed Tomography (CT). The RASER images were acquired in a similar way (Fig. 2D ), but in contrast to the spin echo sequence no RF-pulses were applied. The signal is acquired in the presence of G x and G z field gradients during spontaneous RASER emission, which begins spontaneously, shortly after the crusher field gradient is turned off. For both imaging sequences a crusher gradient is applied to destroy all coherence while negative proton polarization is built up by SABRE pumping at magnetic fields B 0 of 3.9 and 7.8 mT. p-H 2 bubbling is interrupted to allow the solution to settle for a time Δt. For SEI, the image is encoded in the echo signal. In the case of RASER MRI the signal builds up spontaneously in the absence of any RF-excitation. Frequency encoding is performed in x and z direction. 7 The RASER action can be measured over an indefinite period (Fig.3(A) ), when p-H 2 is continuously bubbled through the solution. However, a challenge for imaging under bubbling conditions in the presence of field gradient results from sample motion induced by the bubbling. This collapses the RASER spectrum in each chamber into one average frequency ( Fig. 3(B) ). In order to avoid line-collapse induced by sample motion, and to enable imaging, the p-H 2 flow had to be stopped and an additional waiting time Δt was introduced, which allows for the solution to settle and the motions to halt. Now, both SE and 1D RASER signals could be acquired (Figs.3(D,G) ) shortly after the crusher gradient was switched off. The acquired RASER burst in Fig.3 (G) is significantly longer than the corresponding spin echo in (D) acquired at the same gradient strength of G z = 3.84 mG/cm. The spatial resolution limit is given by δz = w/(γ H G z ) in conventional MRI (22). This limit yields δz SEI = 280 μm for the SEI in Fig.3(E) , and as a result the gap and the edges of the sample are not well resolved. However, for RASER 1D projection in Fig.3(H) , the slope at the image boundaries at the gap is more than three times higher. Thus, a spatial resolution of δz RI ≈ 90 μm is estimated. The measured 1D RASER image in Fig. 3 (H) shows signal lobes outside the boundaries of z = 4 mm, in accord with the simulation shown in Fig. 1(D) . Such artifacts from 1D RASER MRI are analyzed in section 4(b) of the SI and a potential correction method is proposed. Both a 2D-SEI ( Fig.4(A) ) and a 2D-RASER MRI ( Fig.4(B) ) of the same sample are obtained, extending 1D imaging to 2D imaging by reconstructing from 30 angular directions. The field gradient used for the SEI image was 3.5 times larger than that for RASER MRI in order to obtain comparable resolution. Each individual projection in the SEI has a resolution of 50 μm, only about an order of magnitude higher than modern micro imaging (42-44). The two semicircleshaped halves and the 1-mm gap are clearly visible in both images of Fig.4 . These images also display typical projection reconstruction star artifacts outside of the imaging domain. The 2D-RASER image in Fig.4 (B) shows sharper features, but also exhibits a deformed shape of the sample and its gap, paired with several interfering lines. These lines are probably caused by residual motion of the liquid after turning off the pH 2 -pumping. They can be identified in the individual 1D projections, which are used to reconstruct the 2D RASER image (see SI, Fig.2 (C, D) from different angles by varying G x and G z such that G x 2 +G z 2 = const. In (A), the two capillaries used for p-H 2 supply are visible around x = -1 mm, z = 0.5 mm and x = -1.5 mm, z = -2 mm for each chamber. The RASER image (B) is recorded at a 3.5 times smaller gradient than (A), but both spatial resolutions are similar. The RASER image is affected by interference lines. The origin of these artifacts is discussed in the text and in the SI, section 5.

A stark contrast of RASER MRI to traditional MRI is the dependence of RASER MRI images on the magnitude of the nuclear spin polarization. Figure 5 shows a series of 1D RASER and SEI images of the phantom, acquired with decreasing levels of polarization, i.e. decreasing population inversion d 0 . The polarization was adjusted by implementing an increasing waiting time Δt between the polarization step and acquisition.

For SEI, decreasing polarization entails decreasing SNR for each image in Fig. 5 (A), but the shape of the image in the interval 2s < t < 20 s (about a few T 1 relaxation periods) remains invariant. The spatial resolution for the SEI is determined by the slope on the sample boundaries with δz SEI ≈ 50 μm. This observation is in overall good agreement with the theoretical expectation of δz SEI = w/(γ H G z ) = 55 μm. Although the initial negative polarization (d 0 ) changes by more than a factor 10 within the first 20 s, the shape of the SEI images is invariant. This behavior is because the widths of the underlying PSFs barely deviates from a Lorentzian linewidth and radiation damping effects are insignificant. At longer waiting times (t > 20 s), noise 9 becomes more dominant and the shape deteriorates as more efficient relaxation at the walls decreases the image amplitude at the boundaries of the sample.

In contrast, the RASER image shape in Fig. 5 (B) strongly depends on polarization. We attribute the differences between the two image halves to disparities in the bubbling rates and phantom shapes (details see SI section 4(c)). Fig. 5 (C) shows simulated RASER images for five different initial population inversions d 0 and corresponding profiles ρ d (ν) (see SI, Fig. S10 ) to examine the origin of the RASER image distortions. The experiment at Δt = 8 s matches the simulation with only one peak (width = 0.6 Hz, SI, Fig. S9 ) and for the experiments Δt < 8 s, the simulation qualitatively reflects the amplitude deformations and side lobes seen in the measured images. The ripples in some images in Fig. 5 (B) cannot be simulated assuming a uniform division of the RASER image into N = Δ/δν slices. Motional artifacts and variations of T 1 , T 2 * , and B 1 -field over the image domain may be responsible for the observed ripples. 

The proof of principle experiments provided here and the corresponding nascent theoretical framework motivate several new challenges, and promise an opportunity to overcome current fundamental limitations in MRI, potentially leading towards MRI with better signal to noise ratio and spatial contrast. There is absolutely no background signal from other protons (e.g., water) because there is no RF excitation and the RASER signals solely stems from negatively polarized molecules at low concentration (45). The results of this work pave the way to background-free clinical RASER MR imaging with new and more sensitive contrast mechanism. Moreover, the absence of RF excitation (and by extension virtually zero Specific Absorption Rate (SAR)) and the use of substantially lower field gradients requirements (and by extension significantly reduced concerns over peripheral nerve stimulation) in RASER MRI offers new unprecedented standards of patient safety (46). Furthermore, RASER MRI theory is connected to many seemingly disjunct applications in science and technology. The developed system of differential equation Eqs.(1-4) and its solutions for the RASER MRI model are equivalent to the fundamental equations in many other fields (see SI section 6) with prominent examples of synergetics (15) and non-linear dynamics (14, 34, 47, 48) . TT is a founder, equity holder, and President of Vizma Life Sciences LLC (VLS). VLS is developing products related to the research being reported. The terms of this arrangement have been reviewed and approved by NC State University in accordance with its policy on objectivity in research. MSR is a founder and equity holder of VLS. EYC declares a stake of ownership in XeUS Technologies LTD. All other authors have no competing interest.

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Figs. S1 to S11

SABRE samples were prepared under Schlenk conditions. The samples contained 5 mmol/l [Ir(cod)(IMes)Cl] of the SABRE catalyst precursor (41), and c pyr =100 mmol/l of pyrazine in methanol-d 4 . Pyrazine was chosen because it is associated with a single resonance in the NMR spectrum with ns pyr = 4 chemically and magnetically equivalent protons, ideal for RASER and spin echo imaging experiments. For single chamber experiments, the sensitive volume was filled with 900 μl solution, while for the experiments using two chambers with each 300 μl were filled into each side giving a total sample volume V s = 0.6 ml. A glass capillary (~100 μm outer, 30 μm inner diameter), was introduced into each chamber for parallel p-H 2 supply. During polarization build-up, p-H 2 was bubbled though the solution at ~ 30 sccm and at 2 bar pressure. Parahydrogen was generated using a Bruker p-H 2 generator at 35 K, yielding ~ 94% enriched p-H 2 gas. Typically a negative pyrazine proton polarization of P H ≈ -10 -3 to -10 -2 is achieved in a magnetic field ranging from 3.9 to 7.8 mT. This chosen magnetic fields are close to the field B 0 = 6.5 mT, where the SABRE 1 H polarization for pyridine and similar chemical motives such as pyrazine is maximized (18). With respect to RASER MRI low magnetic fields do offer the additional advantage of lower susceptibility artefacts.

A SABRE induced 1 H polarization of P H = -10 -3 corresponds to a population inversion d 0 = c pyr ·V s ·(-P H )·ns pyr = 0.1 mol/l ·6·10 -4 l·-(-10 -3 ) ·4 = 1.4·10 17 . The total number of 1 H spins in the sample is N s = 1.4·10 20 . Analogue calculations yield the initial conditions for simulations in the main text and the SI. For example in Fig. 5 of the main text, the initial population inversion is assumed between d 0 = 3.6·10 16 and 2·10 17 . The 1 H NMR parameters of pyrazine were measured to be T 2 * = 0.7 s (Lorentzian width w =1/(πT 2 * ) = 0.455 Hz). T 1 values at different positions were measured using the results of the SEI images versus t (see Fig. 5 (A)). We found T 1 = 5.0 s in the bulk. The measurement close to the walls varied around T 1 = 2.5 ± 0.5s. For the simulations we chose a difference in T 1 between the bulk and the walls of 3 s. The sample is located in a cylindrical glass tube (ID = 8 mm), divided by a glass slide (1 mm thickness) for twochamber experiments. The designed phantom is hand-made. The 1mm thick glass sheet is held in place by chemically resistant glue. The liquid sample inside the two chambers is located in the sensitive volume of a cylindrical NMR detection coil (ID = 10 mm, height = 10 mm) which is connected to an external resonator with high quality factor (EHQE, Q ext = 360 @ 166 kHz) for sensitive detection of the NMR or RASER signals (49). The total quality-factor of the combined resonator (external resonator and NMR coil) is Q = 100. The B 1 field profile from the NMR detection coil in the center of the sample is calculated to be about 10% lower compared to the field at the edges of the sample. As the RASER active slices interact through the B 1 field of the coil, the coupling now depend on space, which is not accounted for in the parameter β in Eqs.(1-3). In summary, the dependence of B 1 , T 2 * and T 1 on the location of the sample are major sources for RASER imaging artifacts. Correction algorithms for artifacts are state of the art for high field MRI scanners(50) and could mostly be adapted to the artifacts presented here. The magnetic fields of the low cost MRI system are generated by a set of four handmade shim gradients (G x , G y , G z and G crush ) and by an electromagnet producing a constant field of the range of 0.5 mT < B 0 < 20 mT. For our experiments we chose B 0 = 3.9 mT and 7.8 mT corresponding to 166.6 kHz and 333.3 kHz 1 H resonance frequency. The rf frequency of the spectrometer is chosen such that the off-resonance frequency ν 0 is between 20 Hz and 150 Hz away from the 1 H resonance frequency. The homogeneity of the B 0 field is 1 ppm/cm 3 . The p-H 2 supply in a lowfield electromagnet in conjunction with sensitive EHQE detection avoids the necessity of a shuttling system for rapid transport of the sample into a high-field magnet. The G x and G z gradients were used to obtain projections from 30 different angles (in 6 degree steps). All data was acquired in a single scan. Spin echo images were acquired at an echo time of 1 s. 2D images were obtained after projection reconstruction of the 1D slices using a matlab code, written for this project. The spatial resolution is divided into a resolution along a slice in radial and angular direction. The radial resolution is 50 μm for SEI at 21.6 mG/cm, which corresponds to 160 points along the 8 mm sample diameter. The angular resolution with 30 slices spanning 180° is 6°.

There are frequency shifts due to slow magnetic B 0 field drifts in the order of a few ppm per minute. At 333 kHz (7.8 mT) these drifts on a time scale of 10 minutes were more pronounced compared to 166 kHz (3.9 mT). The reason is thermal instability of the current supply in conjunction with heating of the resistive B 0 field coil. For one 1D RASER image measured at 7.8 mT with a corresponding RASER burst lasting a few seconds a drift of a few ppm per minute means less than 0.1 ppm or 0.03 Hz frequency drift. The image domain is typically chosen between 10-100 Hz (corresponding to about 20-200 slices for SEI), so the drift for a single 1D RASER image is negligible. For a 2D RASER image with a total measuring time of about 30 min for all 30 1D slices, the central frequency between the individual 1D slices could differ by a few Hz. Thus, each 1D image was shifted to yield the same center frequency for all 1D images before projection reconstruction.

The simulations based on the model Eqs.

(1-4) were performed using Mathematica 8. The NDSolve[] routine was used for the numerical evaluation of the variables d μ (t), A μ (t), and ϕ μ (t). The computation time of the system Eqs.(S5-S8) can be quite long depending on the number of modes N. All parameters d μ , A μ and ϕ μ are coupled in between each other in a nonlinear way by the cos-and sin-terms on the right sides of Eqs.(S5-S7). This is the reason for many non-linear phenomena which can arise in RASER MRI model, ranging from phase locking, collapse phenomena, non-linear image distortions, edge effects, multiple-period doubling and chaos. While there are exactly N coupling terms for A μ and ϕ μ in Eqs.(S6, S7), the number of coupling terms for d μ in Eq.(S5) is N(N-1)/2. For larger number of slices, N >100, the system of equations becomes elaborate and a large amount of computation is required. In fact the computation time is roughly proportional to N 3 , so the system Eqs.(S5-S7) is classified as a polynomial problem. A typical numerical evaluation using an I5 processor with 8 GB RAM takes about 60 s for N = 50 and can be many hours to days for N > 100.

For these simulations, initial conditions for d μ (0), A μ (0), and ϕ μ (0) are required. The initial conditions for d μ (0) 

This supplement discusses the theory and the physics of RASER modes in the presence of a magnetic field gradient. It aims to explain many phenomena for RASER MRI in one dimension. For this purpose the supplement is subdivided into six subsections. We first introduce a model to extend the nonlinear multi-mode RASER theory (12, 16) , which describes N nonlinear interacting RASER modes (9, 15, 20) . Section 1 describes the derivation of the equations of motion governing RASER MRI in one dimension. A spatial encoding procedure is added which divides the imaging domain Δ = γ H G z L in N nonlinear interacting slices, where γ H , G z and L are the 1 H gyromagnetic ratio, the magnetic field gradient dB 0 /dz, and the sample extension, respectively. In section 2, we derive the point spread function (PSF) for one single RASER mode or slice (N = 1) excluding T 1 relaxation: The hyperbolic secant. In section 3, T 1 relaxation is included and it is shown that the corresponding PSF as a function of time is an asymmetric hyperbolic like secant. In section 4, RASER MRI simulations for N > 20 interacting slices are discussed, especially for a rectangular profile (4a), a rectangular profile superimposed by a sinusoidal modulation (4b) and for a profile close to our experimental two-chamber setup (4c). In section 5, possible artifacts from the 2D and 1D RASER images from the main text are discussed. Finally, section 6 concludes with a discussion of the relation between RASER MRI and other fields of science and mathematics.

As a starting point, one dimensional imaging is discussed. In this case, only one magnetic field gradient G z = dB 0 /dz is applied and the resulting images are one-dimensional projections. We assume that the sample has been polarized into a state of negative spin polarization. In this manuscript, the sample is pumped by SABRE, but nothing precludes the use of other hyperpolarization techniques. We additionally assume, that there is no additional pumping (the pumping period is over) and there was a small but sufficient waiting time to assure that there is no more movement in the sample (see Fig. 2 , main text). This will be the initial time at t = 0 which is characterized by a total population inversion d 0 . In a onedimensional model, the sample is characterized by its extension L along the z-direction, a center z 0 and two boundaries at positions z 0 + L/2 and z 0 -L/2, as shown in Fig. S1(A) . If the distribution of the population inversion at position z is the one-dimensional density function S d (z), then S d (z)dz is the number of negatively polarized spins in the spatial interval dz at the location z. In the presence of the gradient G z the z-dimension is transformed into the frequency space according to a linear transformation ν = γ H ·G z ·z. After this linear transformation the density S d (z) transforms into the one dimensional density ρ d (ν), where ρ d (ν)dν is the number of negatively polarized spins at the frequency interval [ν,ν+dν]. The function ρ d (ν) is shown on the right side of Fig. S1 with a center frequency ν 0 and the two boundaries at ν 0 +Δ/2 and ν 0 -Δ/2. The image domain in frequency space is given by Δ = γ H G z L. Given the two density functions S d (z) and ρ d (ν) in the spatial and frequency domain the total initial population inversion d 0 of the sample is given by the integral

An assumption made for the simulations RASER MRI is that RASER activity in the presence of G z is self-organized. Thus, N >>1 modes or slices are introduced to enable numerical simulations of the theory. A good estimate for the number of slices is given by N = Δ/δν, where the slice separation is chosen δν < w. We found by numerical evaluation that if δν < w, the properties and the shape of the resulting RASER images do not change (invariance principle I). For a simulation of a RASER image without numerical artifacts in the image domain Δ = γ H G z L, a division into N =Δ/δν ~Δ/w as slices is sufficient, where w as ~1/T 1 is the width of the asymmetric point spread function of a single RASER mode (see section 3). Each RASER-active slice with label μ = 1..N can be attributed to a center frequency ν μ given by

The ensemble RASER image can be described by a superposition of N = Δ/δν slices. In contrast to regular MRI with its Lorentzian shaped PSFs, the spectrum associated to each slice is very complicate, because all slices are all coupled between each other. This is shown in section 4(a).

In the following we demonstrate the essential steps leading to the RASER MRI equations, which describe the dynamics of all μ = 1..N RASER slices. We start with the multi-mode RASER theory as described previously (12, 16) . These are given by a set of 2 N non-linear coupled differential equations for N modes. The dynamics of each mode or of the spin species with angular frequency ω μ is described by the population inversion d μ and the complex valued transverse spin component α μ =A μ exp(iϕ μ ), where A μ and ϕ μ are the amplitude and phase of each mode, respectively (12):

The form of these equations, which is formulated in the rotating frame, are based on the multi-mode Laser equations, as published by H. Haken (15) , and by applying the adiabatic elimination of fast variables (the slavery principle) to the fast decaying electromagnetic field modes of the LC resonator (12) . For one single mode (μ = 1) the form of Eqs.(S3,S4) is identical to the extended Bloch equations for the longitudinal and transverse magnetization M z and M T , respectively used in many radiation damping studies (4, 24, 29, 30, 36, 37, 40) The coupling parameter in Eqs.(S3,S4) is given by β = μ 0 ħγ H 2 Q/(4V s ), where μ 0 , ℏ and γ H denote the vacuum permeability, Planck's constant and the 1 H gyromagnetic ratio. V s is the sample volume and Q the quality factor of the resonator, respectively. The coupling constant β is related to the resonator damping rate κ m =ω 0 /Q and to the magnetic coupling constant |g m | 2 = γ H 2 μ 0 ħω 0 /(4V s ) by β =|g m | 2 /κ m . The factors 1/T 1 and 1/T 2 * represent the longitudinal and effective transverse relaxations rates. Each population inversion d μ of mode μ in Eq.(3) is pumped by the rate Γ μ towards the equilibrium population inversion d μ,0 , decays with the rate 1/T 1 and is diminished by the last term on the right side of Eq.(S3), which depends on the sum over all quadratic terms α ν α σ * +α ν * α σ . The transverse modes α μ in Eq.(S4) decay with the rate 1/T 2 * and oscillate with the off resonance angular frequency ω μ . For NMR spectroscopy the origin of these frequencies could be different chemical shifts or splittings due to J-coupling. The last term in Eq.(S4), proportional to the term β d μ α τ , is the source for the RASER emission of mode μ and is responsible for collective phenomena. The term β d μ α μ for τ = μ is responsible for RASER emission. The product β d μ = μ 0 ħγ H 2 Q d μ /(4V s ) can be described as radiation de-damping or damping terms of mode μ, depending on whether the sign of d μ is positive or negative, respectively. This is analogue to radiation damping described in the modified Bloch equations (4, 24, 29, 30, 36, 37, 40) . The cross terms, β d μ α τ for τ ≠ μ, are not included in the extended Bloch equations, but are essential to describe the dynamics of RASER MRI. Their sum determines the time evolution of each of the amplitudes A μ and phases ϕ μ . This has various consequences for RASER images, such as edge artefacts, and non-linear amplitude deformations, as will be shown later.

The RASER MRI equations for one 1 H spin species can be derived based on Eqs.(S3,S4): First, the spin system is pumped into a state with highly negative 1 H polarization P H , i.e. into a state with a large total population inversion d 0 (see Eq.(S1)). For 1 H SABRE pumped organic molecules such as pyrazine or pyridine, the proton polarization is typically in the range from P H ~ -10 -3 up to -10 -1 , which for our experimental conditions corresponds to values d 0 ~10 16 -10 18 (for calculation see Materials and Methods section). To prevent any RASER emission, a strong crusher gradient or a detuned resonance LC circuit ensure that α μ = 0 during the pumping. Thus, only the pumping and T 1 relaxation terms on the right side of Eq.(S3) are relevant during the pumping period. Second, after the pumping and the strong crusher gradient is switched off (all Γ μ = 0) and provided the weak gradient G z for imaging is not too strong, some slices start to be RASER active and oscillate in the presence of G z . The angular frequency ω μ = 2πν μ of each RASER active slice is given by Eq.(S2). After splitting Eq.(S4) into a real part describing the amplitudes A μ and an imaginary part, which describes the phases ϕ μ , a set of 3N nonlinear coupled differential equations for the variables d μ , A μ and ϕ μ is obtained, i.e. The density ρ d (ν) can be seen as a given input profile after pumping of the sample, which depends on the shape of the imaged object. Given ρ d (ν), the quantities d μ , A μ and ϕ μ can be evaluated by numerical evaluation of Eqs.(S5-S7) with the given boundary conditions Eq.(S8), and the measurable total transverse RASER signal results as a superposition of all functions α μ = A μ exp(iϕ μ ).

The Fourier transformation of the total signal Sig(t) results in the RASER image. Note that no boundary conditions for the initial amplitudes A μ (0) and phases ϕ μ (0) are given in Eq.(S8). The reason is that the resulting RASER images are practically independent to the initial values of A μ (0) and ϕ μ (0), no matter whether A μ (0) and ϕ μ (0) are random or defined values. Provided the condition T 1 >> T 2 * holds (which is mostly the case), extensive numerical simulations reveal three invariance principles with respect to the Fourier transformed RASER images in the absolute mode: I. The RASER image does not depend on to the slicing δν as long as δν < w. II. The RASER image contrast and resolution is independent of T 1 , and III. The RASER image is invariant with respect to the initial conditions A μ (0) and ϕ μ (0). These three invariance principles have significant consequences for RASER MRI. The invariance from the choice of the slicing distance δν (principle I) means that a continuous limit N   exists, where the discrete variables d μ (t), A μ (t) and ϕ μ (t) become continuous variables. In this limit all sums in Eqs.(S5-S8) become integrals. Another consequence is that numerical simulations produce reliable results without changing the involved physics as long as δν < w. We found that ripple like artifacts arise in the image if δν ~ w but the envelope of the image is the same compared to δν << w. However, if δν << w, the numerical simulations can become very time consuming. A good compromise is δν ≈ 1/T 1 ≈ w as , where no ripples are visible. The invariance of the RASER image shape and contrast on the value of T 1 , (principle II), means that the contrast of RASER MRI cannot be associated to the width of single PSFs. This includes the width of the asymmetric PSF w as ≈ 1/T 1 of a single RASER mode, introduced in section 3. The contrast mechanism is based on collective interactions, as will be shown in chapter 4(a). The invariance of the RASER image on the initial conditions A μ (0) and ϕ μ (0) (principle III) allows for reproducible results irrespective of the noise excitation. The amplitude and shape of the RASER image do not change if the RASER burst is initiated either by spin noise or a defined excitation sequence with a weak RF or DC pulse. We found that for very different initial conditions A μ (0), ϕ μ (0) (more than one order of magnitude), there is a small shift of the entire RASER burst signal. This leaves the absolute Fourier transformed spectrum invariant but the phased spectrum is shifted by a global phase. This small shift can be avoided by applying a weak DC or rf pulse to initiate the reproducible RASER bursts, beneficial for averaging and 2D RASER MRI.

A detailed mathematical analysis of the dynamics of the image formation and how to explain the three invariance principles is quite elaborate and not the main focus of this contribution. Therefore, we focus on numerical simulations to investigate image artifacts and the more sensitive contrast mechanism. In section 4(a) we use the simple example of a rectangular profile. Due to the invariance principle III, we chose the simplest case, ϕ μ (0) = 0 and a constant small value A μ (0) ~ 10 9 -10 10 . The value for A μ (0) is in the order of the spin noise amplitude N s 1/2 (26), where N s ~10 19 -10 20 is the total number of spins in the sample.

Finally, we analyze whether a local threshold condition exists for RASER MRI. Without gradient, there is only one RASER mode (N = 1). The threshold condition for RASER action is 0 ≥ th = 4 / ( 0 ℏ H 2 2 * ), (12) where d th is the threshold population inversion. An equivalent expression for this threshold condition is ε = d 0 /d th = T 2 * /τ rd  1, where ε is a dimensionless quantity characterizing the 21 threshold. However if a gradient is applied, the sample divides into N = Δ/δν >> 1 slices and the threshold condition ε = d 0 /d th > 1 cannot be used any more as a proper criterion for RASER activity.

When a gradient is applied, the population inversion d 0 is distributed over the image domain Δ according to the population inversion density ρ d (ν). Assuming no interaction between the slices δν, a threshold condition can be formulated based on the population inversion within a frequency interval w = 1/(πT 2 * ). Within this frequency interval the damping with rate 1/T 2 * is compared to the RASER dedamping process. A region r μ = [ν μ -w/2, ν μ +w/2] at position ν μ and width w = 1/(πT 2 * ) is RASER active if be RASER active while the regions at the boundaries are not. Close to the center the slices cooperate with all its neighbors in a constructive way, while the slices close to the two boundaries cooperate destructively with their neighbors. Because of the cooperative action between slice μ with all other slices, a local threshold condition does not exist and has to be replaced by a non-local threshold condition. Numerical simulations show that RASER action for a region r μ in the domain Δ depends on d 0 , on the width of the image domain Δ and on the detailed shape of ρ d (ν). A detailed mathematical evaluation of this statement is quite elaborate, but a numerical evaluation of the case using a rectangular profile is shown in section 4(a). In nearly all simulations shown in the following the threshold population inversion density 

We will proceed with the simplest possible case, and derive the exact solution for one RASER mode (N = 1) if T 1 relaxation can be completely neglected (T 1 = ∞). In this case an exact solution in the form of tanhand sech-functions exists. The relevance for the line shape in NMR spectroscopy due to strong radiation damping effects has been demonstrated (4, 24, 29, 30, 35-37, 40) . Here, we repeat the derivation of the basic results described there for two reasons. The first reason is to be compatible with the nomenclature used here. Secondly, Mao et al. (29, 35, 36) studied amongst others the radiation damped signal burst after a radio frequency pulse excitation angle close to π, which is in close correspondence to one self induced RASER burst in Eqs.(S5-S8). For N = 1 and neglecting T 1 relaxation Eqs.(S5-S8) are reduced to its most basic form, 2 * 2 1 4 , (S10), 1 , (S11), , (S12).

We choose as boundary conditions for Eqs.(S10-S12)    

Eqs.(S11,S12) are equivalent to one single differential equation for the complex valued transverse spin component

. The derivation of the exact solution for d(t) and A(t) follows from differentiation of Eq.(S10), which results in 22

After substituting the right side of Eq.(S11) for At  in the preceding expression we obtain 2 2 * * Eq.(S10) 22

The solution is given by integration, which results in

The Integration constant C is fixed by the initial conditions A(0) ~ 0 << d 0 , d(0) = d 0 , so dd/dt = -4βA 2 (0) ~ 0, giving C = d 0 (2/T 2 * -β d 0 ). After introducing the parameter q = (1+ ε 2 -2ε) 1/2 , where the dimensionless variable is ε = T 2 * /τ rd = T 2 * β d 0 , the previous equation for U can be written as 22 ** 22

Introducing the parameters a = 1/(T 2 * β) and b = q/(T 2 * β), Eq.(S14) simplifies to

which is equivalent to

This integral is known as 

The exact form for the transverse spin component A(t) can be derived by differentiation of Eq.(S16). Using the relation

Eq.(S17) can also be written as

Applying the square root on both sides of the preceding equation, A is finally

According to Eq.(S18), A is a hyperbolic secant function (soliton solution in the time domain) and is symmetric with respect to t 0 , which is the time of maximum amplitude. An interesting feature of Eq.(S18) is that the envelope of the Fourier transformed spectrum in the frequency domain ω is once again a sechfunction. According to Mao et al. βd 0 ≈ 1 the factor q << 1 becomes very small, so the associated linewidth w sech of the corresponding spectrum Eq.(S22) is much smaller compared to the linewidth w = 1/(πT 2 * ) of a Lorentzshaped peak, the latter representing the standard PSF for Spin Echo Imaging (SEI). The full width at half maximum w sech is determined by the argument [πT 2 * ω /2q] in Eq.(S19), which for ε ≈ 1,  0 ≈ π and w = 1/(πT 2 * ) is w sech = (2/π) ln(2+3 1/2 )(ε -1) w = 0.84(ε -1)w. For ε > 2.2 the width w sech > w, while in the range 1< ε < 2.2 closer to threshold w sech < w. For example, at ε =1.4 w sech = 0.336 w = 0.15 Hz for T 2 * = 0.7 s. The argument w sech < w for 1< ε < 2.2 holds even for the case of finite T 1 relaxation associated with an asymmetric PSF, as will be shown in the next section. Fig. S2(A) . The four point spread functions (PSF`s) are displaced in the frequency dimension (with center frequencies at 48, 54, 57 and 58 Hz) for clarity. Note that in each PSF and with decreasing ε the spectral width and the spacing between subsequent oscillations is decreasing. Fig. S2(A) shows the numerical simulations of four different PSF for the case T 1 = . The transverse spin component A(t) versus time t is plotted for four different parameters ε = T 2 * /τ rd = T 2 * β d 0 , ranging from ε = 7 (green) far above threshold, ε =2.8 (red), ε =1.55 (blue) to ε =1.16 (orange) close to threshold. The simulation parameters are T 2 * = 0.7s, Q = 100, V s = 0.5 cm 3 and the initial conditions are ϕ(0) = 0, A(0) = 10 13 . All four simulated RASER bursts are symmetric sech-functions which are identical to the analytical solution given by Eq.(S19). Furthermore, each PSF is symmetric with respect to the time of maximum signal, t 0 . The time t 0 , as given by Eq.(S20), and the width of each burst increases for decreasing values ε. As ε approaches the threshold (ε = 1) the duration or width in the time domain of each PSF becomes 25 arbitrarily long, and the maximum amplitude arbitrarily small, but the area below each of the RASER envelopes remains the same. The corresponding Fourier-transformed spectra are shown in Fig. S2(B) . The four PSF`s are displaced in the frequency space for a better separation. Note that with decreasing ε the spectral width and the period due to the modulation proportional to the factor cos(ωt 0 ) (see Eq.(S19)) in each PSF is decreasing. Close to threshold the corresponding PSF are much narrower compared to the linewidth of a conventional Lorentz-peak. In our experiments typical measured values are T 1 ~3 s -5 s and T 2 * = 0.7 s, so the PSF has a width in the range 0.2 Hz < w as < 0.33 Hz and the Lorentzian line width is w = 1/(πT 2 * ) = 0.455 Hz.

Case N = 1 including T 1 relaxation: The asymmetric Point Spread (a-PSF).

For this next case, we include T 1 relaxation. To our knowledge this has not been discussed in the literature. The dissipation caused by T 1 relaxation, represented by the additional term -d/T 1 , is introduced in Eq.(S10) for one mode (N =1). No analytical solution exists for this case, and the symmetry of the PSF given by Eq.(S18) with respect to time t 0 is lost. Fortunately the key properties of the corresponding asymmetric PSF can be evaluated by numerical simulations of Eqs.(S10-S12), including the additional loss term -d /T 1 . Fig. S3(a) . Far above threshold, at ε = 7, the peak amplitude including the modulation due to the cos(ω t 0 )-term of the spectrum looks quite similar to the spectrum for T 1 = ∞ in Fig. S2(b) . Close to threshold, at ε = 1.4 and 1.08, the peak amplitude is tens of times smaller compared to ε = 7 and there is nearly no phase modulation visible.

To find out how the a-PSF depends on the parameter ε we performed various simulations as shown in Fig.  S3 , where a series of six PSF in the time and frequency domain are shown for different values of ε. In Fig.  S3 (A) the envelope of the amplitude of A is plotted versus time for ε =1.08 (black), 1.4 (blue), 2 (brown), 2.8 (red), 4.4 (orange) and 7 (green). The simulation parameters are given in the caption of Fig. S3 . Due to T 1 relaxation, all RASER bursts look like slightly asymmetric sech-functions (a-PSF). The weak asymmetry is with respect to their maximum signal at time t 0 , i.e. the duration of the initial rise of all signals is slightly shorter compared to the decaying tail. T 1 relaxation is also responsible for large differences in the maximum amplitudes between each a-PSF. Close to threshold at ε =1.08, the maximum amplitude is several hundred times smaller in comparison to ε = 7, far above threshold. Conversely, the duration of the a-PSF at ε =1.08 is tens of times longer when compared to the duration at ε = 7. Fig. S3(B) shows a plot of the corresponding phased Fourier transformed spectra of the six a-PSF in Fig. S3(A) .

Far away from threshold, at ε = 7, the peak amplitude and the modulation due to the cos(ωt 0 )-term of Eq.(S22) of the spectrum looks quite similar to the spectrum for T 1 = ∞ in Fig. S2(B) . At ε = 1.4 and 1.08 close to threshold, the peak amplitude is tens of times smaller compared to the peak amplitude at ε = 7. Additionally, close to threshold there is nearly no phase modulation visible in the spectrum, in contrast to the phased spectra close to threshold in Fig. S2 (B), which are characterized by a narrow envelope modulated by several oscillations. This difference in the number of periods visible in the PSF in Fig.  S2 (B) and a-PSF in Fig. S3 (B) close to threshold is directly connected to the time of maximum amplitude, t 0 . To analyze this key property, the time t 0 as a function of ε has been numerically evaluated for different values of T 1 . The result is shown in Fig. S4 , where the simulation parameters are given in the caption. For the case T 1 = ∞, t 0 decreases monotonically with increasing ε. Note the dotted curve, which represents the exact solution given by Eq.(S20), is in good agreement with the numerical solution given by the solid blue line. A different behavior is observed for finite values of T 1 , here ranging from T 1 = 20 s, 10 s, 6 s to 3 s. Starting at ε = 1, t 0 increases with increasing ε until a maximum value is reached in the range 2 < ε < 4 and finally t 0 (ε) decreases until approaching the curve for T 1 = ∞ in an asymptotic way. Here we state the similarity between the PSF and the a-PSF far from threshold and significant differences close to threshold. 27 

Until now, we have studied PSFs for one single slice (N = 1). To describe RASER MRI, the concept of a single PSF is insufficient. Thus, in the following subsections, we will explore the physics of 1D RASER images consisting of many interacting slices (N ≫1). Now collective effects dominate the physics of image formation and are essential to describe the image contrast as well as the nonlinear artifacts that arise. To understand these phenomena, we chose simple profiles to reproduce typical RASER MRI features in simulations (see subsections 4 (a-c)).

In all the simulations shown here, we assume parameters matching our experimental conditions: T 1 = 5 s, T 2 * = 0.7 s (line width w =1/(πT 2 * ) = 0.455 Hz), quality factor Q = 100, a cylindrical sample with volume V s = 0.5 cm 3 and diameter L = 0.8 cm. All these parameters have either been measured in our 1 H RASERimaging experiments at 166 or at 333 kHz 1 H Larmor frequency (B 0 = 3.9 or 7.8 mT) with SABRE pumped pyrazine, or correspond to parameters of the setup. The typical range for the magnetic field gradient is 210 -4 G/cm < G z < 210 -2 G/cm, which for a sample dimension of L = 0.8 cm corresponds to image domains in frequency space ranging from 0.67 Hz < Δ < 67 Hz. Three different spin density profiles are defined in Fig. S5 : A rectangular profile above threshold, a profile with a constant value superimposed by a sinusoidal profile slightly below threshold, and a profile which is close to the experimental conditions with two sections with flat tops separated by a gap and with edges described by tanh functions. (c) Two regions or sections separated by a 1 mm broad gap, to mimic the experimental conditions. The profile is approximated by tanh-functions, i.e. ρ d (ν) = 4.5•10 15 /Hz (0.5(tanh[(ν-ν 0 + Δ/3)/0.7] -tanh[(ν-ν 0 +Δ/15)/0.1]) + 0.4 (tanh[(ν-ν 0 -Δ/15)/0.1] -tanh[(ν -ν 0 -Δ/3)/0.7])) (blue). With a chosen value for the image domain Δ = 10 Hz and w as = 0.33 Hz this corresponds to N = Δ/w as = 30 slices. The dotted line at ρ d th = 6.6•10 15 /Hz indicates the population inversion density at the RASER threshold, which is related to the threshold population inversion d th = ρ d th w = 3•10 15 . Note that for the rectangular profile in Fig. S5 all values fulfil ρ d (ν) > ρ d th , while for the cosine-modulated profile ρ d (ν) < ρ d th . As will be shown in section 4(a) this does not mean that no RASER activity is observed for the cosine-modulated profile.

The simulation of a RASER image based on a rectangular profile with constant ρ d (ν) = ρ d rect has two different purposes. First, it serves as a simple model system to analyze nonlinear phenomena. Second, it serves as a reference to correct for nonlinear amplitude deformations, which occur in arbitrarily shaped RASER images.

An overall view of how the amplitude and shape of a RASER image depends on different values of the initial population inversion d 0 is shown in Fig. S6 . In S6(A), panels I-V, five simulated RASER signals are shown for five different values d 0 ={5.5, 6.6, 10, 15, 30}10 16 . As before, the rectangular profile is centered at ν 0 = 50 Hz and with a constant population inversion density ρ d (ν) = d 0 /Δ in the frequency range 45 Hz < ν < 55 Hz (Δ =10 Hz). With a given step size δν = w as = 0.2 Hz this results in N = 50 slices. The green dotted line in Fig. S6(B) indicates the threshold population density ρ d th =d 0 /Δ = 6.610 1 /Hz. This means that the population inversion within a slice of a width w (ρ d (ν)w) is smaller than the threshold population inversion d th required to start RASER activity. Therefore, at first glance no RASER action is expected for all N=50 slices. Surprisingly, a RASER burst is visible (burst I in S6(A)) and the corresponding image I in S6(C)) is a symmetrical peak with a maximum at 50 Hz and a width much narrower than Δ. The reason for this shape is the cooperative action between all 50 slices. The slices close to the center at 50 Hz successfully surpass the threshold through the cooperation with their neighbors. In contrast to that, the slices on the edges (at 45 and 55 Hz) do not surpass the threshold and are not RASER active. Therefore the simulated image amplitude is maximal in the center of the image domain.

The next case for d 0 = 6.610 16 is shown in image II in S6(C), where ρ d (ν) = ρ d th holds. Now, all slices in the image domain Δ are RASER active and contribute to the image. Nonetheless, the image does not have a rectangular shape, but a bell shaped image with a maximum amplitude in the center is formed. Furthermore decaying sidelobes left and right from the image boundaries arise. We found that cooperative action between all slices leads to broad and complicate spectra of each slice and are responsible for the signal outside the image domain Δ (see Fig. S7 ).

In image III (d 0 = 110 17 ) the amplitude at the edges is nearly as high as the amplitude in the center and the decaying sidelobes are more pronounced. In this case, the slices at the edges of the image domain benefit from the enhanced cooperative interaction between all slices. This effect is even more pronounced in the images IV (d 0 = 1.510 17 ) and V (d 0 = 310 17 ), with ρ d (ν) >> ρ d th , where the amplitude in the center is comparable or smaller in size with respect to the amplitude at the image boundaries. The exact quantitative description and how a rectangular RASER image is formed is quite involved and will be the subject of an independent publication.

Next, we study the sensitivity of a RASER image with respect to small disturbances in the polarization distribution, in the profile ρ d (ν). Therefore we compare simulated RASER images with standard imaging based on Lorentz shaped PSFs. As input profile, we assume a rectangular input profile ρ d (ν) for d 0 = 910 16 and for an image domain ranging from 45 Hz to 55 Hz (Δ = 10 Hz, δν = 0.151 Hz, N = 66). On this rectangular profile, a small perturbation is added in the form of a rectangular shaped hole in the center ν 0 = 50 Hz. The width of this hole is one half of the natural linewidth, i.e. w/2 = 1/(2πT 2 * ) = 0.23 Hz, and the amplitude is 20% smaller compared to the maximum of ρ d (ν). In this case, the smallest value in the hole is still above the threshold population density ρ d th = 6.610 15 /Hz and therefore all N = 66 slices are RASER active. The resulting simulated RASER image is depicted in Fig. S7 (B) as I (black). Comparing both images in Fig. S7(B) , there are two apparent differences in the RASER image: First, similar artifacts as shown in panel III of Fig. S6(C) for a similar population inversion arise. The amplitude of the RASER image with respect to ρ d (ν) inside the imaging domain is deformed and sidelobes outside of the image domain are more pronounced. Both these artifacts are typical for a rectangular input profile above the threshold density and are not attributed to the introduced small perturbation. Second, there is a minimum in intensity in the center of the image at ν 0 = 50 Hz. In the RASER image, the amplitude of this perturbation is about three times lower and the slope at the edges three times steeper. This simple example demonstrates the advantage and problems of RASER imaging, which is sensitive to small perturbations, but suffers from amplitude deformations inside of Δ and side lobes outside of Δ.

Spectra of seven representative slices are chosen for Fig. S7(C) . Each of these spectra are obtained after Fourier transform of the respective RASER signal of slice μ. The slices at the image boundaries (μ = 1 and 66) are depicted in brown, slice μ = 33 close to the center is depicted in blue, while the slices between these extremes are shown in orange (μ = 12 and 54) and green (μ = 24 and 42). The peak amplitude of the spectrum in the center (μ = 33, blue) is significantly smaller compared to slices μ = 24 and 42 (green). This motivates the increased sensitivity of RASER imaging with respect to small changes in the amplitude of ρ d (ν). The width of each spectrum in Fig. S7 (C) is close to w. Furthermore all spectra feature broad sidelobes left and right from the peak maxima, which can extend out of the image domain Δ. The highly resolved local hole with w/2 width cannot be explained by the width w for the individual N = 66 slices. Instead, the collective interaction between all slices generates the local image contrast. This interaction involves collective phenomena such as synchronism, line collapse and other non-linear phenomena outside of the scope of this manuscript, which focuses on the main aspects of RASER MRI and aims to keep the model as simple as possible. 

In the next step we will study simulations of RASER images with more complicated profiles. One interesting example is shown in Fig. S8 , depicting simulations of the time dependent RASER signals 31 (A,D,G) and the corresponding RASER images (B,E,H) for a given offset superimposed by a cosinemodulated profile (red lines in insets of (A,D,G)). The simulation parameters and boundary conditions are given in section 5 and in the caption of Fig. S8 . The threshold population density, ρ d th = 6.6•10 15 /Hz is indicated as green dashed line in the insets of Fig.  S8(A,D,G) . In order to resolve all the details, N = 50 slices are assumed for the simulations, which corresponds to an image domain of Δ = N • w as = 10 Hz. The four periods of the profile have a modulation depth of 20%, i.e. ρ d (ν) = C 0 (1+ 0.2 cos[8π(ν-ν 0 )/Δ]), where the factor C 0 =d 0 /Δ quantifies the average spin density. Three different offset values of C 0 are chosen, as indicated in the insets (A,D,G). In the first case, Fig. S8(A) , C 0 = 8.6•10 15 /Hz, so ρ d (ν)  ρ d th holds inside the whole image domain. The Fouriertransformed image in the absolute mode in S8(B), red line, reflects the four periods of the profile, and the rough envelope of the image with a maximum amplitude in the center at ν 0 = 50 Hz is deformed relative to the profile. We know already from the previous section, that the image of a rectangular profile is deformed and characterized by a Gaussian like shape, reaching maximum amplitude at the center.

The 3He nuclear Zeeman maser

Spin-exchange-pumped 3He and 129Xe Zeeman masers

3He maser for earth magnetic field measurement

Observation of noise-triggered chaotic emissions in an NMR-maser

Solid-State Nuclear Spin-Flip Maser Pumped by Dynamic Nuclear Polarization

Generation of high-frequency current by the products of a photochemical reaction

Spontaneous emission of NMR signals in hyperpolarized proton spin systems

A DNP-hyperpolarized solid-state water NMR MASER: observation and qualitative analysis

A Magic Angle Spinning Activated 17O DNP Raser

Deviation from Berry's adiabatic geometric phase in a Xe131 nuclear gyroscope

Nuclear spin gyroscope based on an atomic comagnetometer

From LASER physics to the parahydrogen pumped RASER

Nuclear Spin Gyroscope based on the Nitrogen Vacancy Center in Diamond

Nonlinear Dynamics and Chaos: With Applications to Physics

Synergertics: An Introduction

SABRE and PHIP pumped RASER and the route to chaos

Transformation of symmetrization order to nuclear-spin magnetization by chemical reaction and nuclear magnetic resonance

SL greatly acknowledges the RWTH Aachen University for accepting me as a guest scientist, providing the research environment and equipment to run all experiments of this study at the ITMC (Institute of Technical and Macromolecular Chemistry). SA thanks his wife Jenny Oquendo Mora for keeping up the moral during difficult situations in the COVID 19 pandemic. Excellent cooperative and IT support from Stefan van Waasen, Michael Schiek and Ulrich Probst from Forschungszentrum Jülich is greatly acknowledged. Michael Adams is greatly acknowledged for valuable help in designing the phantom. 

A correction of the deformed image of the cosine-modulated profile is possible by using the corresponding image of a rectangular profile. We call this procedure a global rectangular correction. The rectangular profiles are drawn as black lines in the insets of Figs. S8(A,D,G), and the constant value corresponds to the average value of the cosine-modulated profile, i.e. ρ d rect (ν) = C 0 . The black lines in (B,E,H) represent the reference images based on the corresponding three rectangular profiles in the insets of (A,D,G). The global rectangular correction procedure consists of dividing the amplitude of the cosinemodulated image by the amplitude of the reference image for each frequency. The results of this procedure are shown as blue lines in Figs. S8(C,F,I). The images of the cosine-modulated profile corrected in this manner can be directly compared to the normalized profiles (red lines in (C,F,I)), the latter being defined by ρ d * (ν) = ρ d (ν)/C 0 . In Fig. S8 (C) there is good agreement between the corrected image (blue line) and the normalized profile (red line), in the image domain 45 Hz <ν < 55 Hz.A case where half of the population inversion densities ρ d (ν) are below the threshold density ρ d th is shown in Figs One might think that RASER activity is impossible. This would be a misconception however, due to the fact that nonlinear interaction can increase the population inversion of those slices near the maxima of ρ d (ν) close to threshold. The corrected image in Fig.S8 (I) retains roughly the shape of the normalized profile, but at the cost of a large deformation in the amplitude, with the minima of the corrected image being much closer to the normalized profile than the maxima.We found that in general for small modulation depths or variations close above threshold, the images reflect the normalized profile. At higher modulation depths and further away from threshold, the nonlinear amplitude deformations increase significantly. For example, for a cosine-modulated profile with 50% modulation depth and ρ d (ν) < ρ d th , the image transforms into one dominant peak at the center and very small peaks close to the local maxima of ρ d (ν). One open question for RASER MRI applications is whether there is an algorithm capable to correct for all nonlinear deformations.

In the last example we attempt to come closer to our real experiment, in which two sections of a cylindrical sample with 8 mm inner diameter are separated by a gap of 1 mm size and pumped separately by SABRE. We determined the 1D projection experimentally with hyperpolarized high resolution Spin-Echo Imaging (see Fig.5 (A), main text). These 1D images can be approximated using two step-like functions separated by the gap and the rising and falling edges of each step by tanh-functions. Examples of normalized two chamber profiles are sketched in the insets of Figs. S10(A). All five profiles correspond to an analytical expression given by the sum over four tanh step functions, i.e. ρ dThe constants A 1 and A 2 define the maximum amplitudes of both steps, ν 0 is the frequency offset, b 1 , b 3 (b 2 , b 4 ) denote the two positions of the rising (falling) edges of the tanh step functions relative to ν 0 , and w 1 , w 3 (w 2 , w 4 ) define the widths associated to the rising (falling) edges. The normalization of ρ d * (ν) to the value one is hereby related to the larger value of A 1 or A 2 , respectively. The amplitudes of A 1 or A 2 are not necessarily equal, but depend on the pumping conditions and T 1 relaxation rate of each of the two chambers. The maximum value of the profile ρ d s. The Fourier transformed spectrum in the phased mode (S9(B)) has the shape of a phase modulated sech function, as discussed in section 3. The corresponding absolute phased spectrum in S9(C) is a symmetrical peak centered at 127.7 Hz and the width at half height is 0.6 Hz. This width is broader compared to the Lorentzian linewidth w = 0.45 Hz. Taking into account the result for the line shape and width for single slices in Fig. S7 (C) this indicates that only a few cooperating coupled slices are responsible for the shape and width of the observed spectrum. The simulated RASER burst based on the chosen ρ d (ν) is nearly symmetrical with respect to its maximum amplitude at t max = 1.5 s (see S8(D)). Both this signal and the phased and absolute spectra in S8(E,F) are in good agreement with their experimental counterpart. Note that the absolute spectrum in (ν) in the insets of S10(A) differ in shape and amplitude ratio between the left and right half of the phantom, which takes into account two different overall T 1 relaxation rates (T 1 = 5 s, 3 s for the left and right half, respectively) as well as locally enhanced relaxation rates at the walls of the sample chambers. The five RASER spectra I-V in Fig. S10(B) represent the corresponding images, which are obtained from the RASER signals in S10(A) after Fourier transformation in the absolute mode. With increasing d 0 , the images change from (I) a narrow peak, to (II) a broader image centered at 127.7 Hz, (III) an image where the right half starts to be visible, (IV) an image where both halves with equal amplitudes separated by the 35 gap appear and finally (V) to an image where the right half is larger in amplitude. The simulated image I in S10(B) is identical to Fig. S9 (F) and is discussed above. The other four simulated images II-V roughly reflect the features of the experimentally measured RASER images in Fig. 5B , main text, including the side lobes in images IV,V appearing outside the image boundaries and non-zero values in the gap.There are two additional phenomena, sometimes observed in the measured images: Pronounced ripples and regions with strongly changing amplitudes. For example ripples can be seen in Fig. 5(B) , main text, at Δt = 1 s, 4 s and 5 s at positions 1 mm < z < 2 mm. An example for peaks with strongly changing amplitudes is shown in Fig. 5(B) , main text, for Δt = 1 s at position z = -1 and 3 mm. Possible imaging artifacts for RASER MRI are discussed in the next section.

This section focuses on the artifacts that can arise in 1D and 2D RASER images. The manuscript mainly discusses 1D images, also called projections, but in Fig. 4(A) and (B) of the main text two different 2D images are shown. A spin echo image for reference and a RASER image. These 2D images are generated using projection reconstruction of 30 1D images measured from 30 angular directions. In these 2D images, artifacts arise due to the projection reconstruction algorithm as well as within the 1D projections themselves.Projection reconstruction generates star artifacts. These are well known, as projection reconstruction is widely used e.g. in imaging methods such as computed tomography (CT). The star artifacts are more pronounced further away from the center. This can for example be seen in the bottom left corner of Fig.  S11 (A) and the top left and right corner of Fig. S11(B) . Star artifacts can be reduced and the angular resolution increased by measuring more projections.Other artifacts and features in the 2D images stem from the p-H 2 delivery system, necessary for SABRE pumping. The p-H 2 is introduced through a capillary in each of the chambers. In the SEI, they can be seen as dark spots each corresponding to the location of a capillary for p-H 2 supply. For RASER MRI, the capillaries can additionally result in nonlinear distortions of a given projection, as RASER MRI is very sensitive to local fluctuations in polarization. Additionally, the p-H 2 bubbling introduces a motion of the liquid in both chambers. This motion stops after about 1-2 s. To ensure a motion-free reference image, Δt = 5s is chosen for SEI. The RASER image however, is very sensitive to the initial population inversion as visualized in Fig. 5 of the main text. To ensure that both chambers have a population inversion in a regime where an image is formed, Δt = 2s is chosen for RASER MRI. This leaves the RASER image more susceptible to small residual motion.Most artifacts can be identified in the 1D projections. Thus, for the 2D RASER image in Fig. 4(B) of the main text, five projections (I-V) are selected. Their gradient directions are drawn as colored lines in Fig.  S11(B) . The 1D images of these chosen angles are depicted in Fig. S11(C) . When choosing a gradient direction perpendicular to the gap between the half circles of the phantom, the gap can be identified as a minimum amplitude in the middle of the projection. Here, projections III (green) and IV(blue) are close to this condition. For the 1D images in Fig. 3 of the main text, the gradient is chosen perpendicular to the gap as discussed there. A gradient parallel to the direction of the gap yields a projection without a minimum amplitude in the middle, similar to the projection of a circle (see Fig. S11 (C), projection V, orange).The most prominent artifacts in the 2D RASER image are interference lines through the entire image. They can be identified within the 1D projections that have a gradient direction perpendicular to the observed interference line. In the projection they can be identified as "spikes" at the given position. This is visualized exemplarily for two interference lines, marked by stars and arrows in the 2D image. They 36 stem from projections II (red) and III (green), respectively. These artifacts can also be identified in the projections and are encircled and marked with stars in Fig. S11(C) . One possible reason for such interference artifacts is the sensitivity of the coupled RASER modes to local disturbances. Disturbances can be caused by residual motion as described above, the capillaries for p-H 2 delivery and fluctuations in the profile ρ d (ν), produced by the SABRE pumping. Further artifacts that arise in a 1D RASER image are leaking signal into a gap as well as sidelobes that arise outside of the image domain. They are small in the 2D images depicted here, but can play a major role in 1D RASER images recorded at other experimental conditions. These leaking and sidelobe artifacts are discussed extensively in section 4a for a rectangular profile ρ d (ν).Other phenomena such a global and relative dipolar shifts might contribute to the observed artefacts. Past studies showed that global dipolar shift effects are quite small for SABRE pumped samples. For example the total 1 H dipolar shift generated by the SABRE pumped liquid, as described in the supplement of (19), was in the order of 200 mHz (B dip ~5 nT). This global dipolar field decays with time as the population inversion d 0 is depleted during the RASER burst. This time dependent global dipolar field may alter the shape of the image by up to 0.2 Hz, which is in the order of the width of one slice δν = w as . Whether these time dependent dipolar shifts contribute to the interference artifacts or could produce chaotic regions is still an open question. To keep the RASER MRI model as simple as possible, we also neglected dipolar contributions in the simulations. 37 

In conclusion we have presented the theory of RASER MRI, which in the imaging domain Δ can be approximated by a network of N = Δ/δν equidistant nonlinear coupled slices (Eqs. (S5-S8) ). Various nonlinear effects, such as high sensitivity to local variations in the input profile ρ d (ν), amplitude deformations and edge artefacts, are predicted by the theory and indeed have been observed in various experiments.Interestingly, the mathematical form of Eqs.(S5-S8) is closely related to many descriptions of nonlinear and collective phenomena in other fields of Science. Especially self-organized processes, the topic of synergetics uses order parameters and adiabatic elimination of fast variables to derive the LASER equations. Looking at the multimode case, H. Haken has shown a close relationship between the LASER assuming a continuous number of modes with the Ginzburg-Landau theory of superconductivity (15) . A similar correspondence between uniformly distributed oscillators with non-local coupling with the Ginzburg-Landau theory of superconductivity is shown in the article by Y. Kuramoto and D. Battogktokh (33). In fact, the Kuramoto model of synchronized oscillators, which is equivalent to Eq.(S7), is a subset of our model of RASER MRI. As outlined in the book from S. Strogatz (14) and in many articles (33, 34, 51) phase synchronism in non-linear coupled oscillators occur in many very different fields in physics and biology. Examples are biological rhythms, synchronized action of fireflies and tree frogs, electrically coupled Josephson arrays (14) , mechanically coupled metronomes (chimera states) (52) spin torque nanooscillators (47) and synchronized spin-valve oscillators (53), and the dynamics of neural oscillator networks (48).Furthermore Hakens equations for one single LASER mode without applying the enslaving principle are apart from a variable transformation identical to the pioneering Lorenz equations (14, 54) , the ladder of which describes chaos in a three dimensional space (15) . In a similar manner it has been shown that chaos arises for the case of two enslaved RASER modes (16) , which involves the evolution of four independent parameters d 1 , α 1 , d 2 , and α 2 . At present the exact analysis for why and how chaos arises in N coupled RASER modes is unknown.Due to the strong links of Eqs.(S5-S8) to many different fields we believe that the presented theory of RASER MRI may be a base for a deeper understanding of self-organizing processes based on both adiabatic elimination of fast variables and on synchronism.