The suggested double angle stimulated echo (DA-STE) sequence including the relevant gradient moments is shown in Fig. 1(a). An extended phase graph [22, 23] diagram is shown in Fig. 1(b) to visualize the spin and stimulated echo formation. Generally, three RF pulses will generate up to four spin echoes and one stimulated echo depending on distance and dephasing [22, 23]. If the refocusing flip angle is not exactly 180°, the magnetization dephased by the first gradient moment will partly be ‘stored’ as modulated longitudinal magnetization [22, 23]. Longitudinal magnetization stays unaffected by gradient moments but relaxes with the time constant T1. The third RF pulse flips the modulated longitudinal magnetization back into the transverse plane to form a stimulated echo [22, 23]. In order to ensure proper separation of the stimulated echo path from the spin echo path, the spin echo is dephased by twice the dephasing gradient moment from the first period [24]. This is emphasized in the extended phase graph schematic in Fig. 1(b).

a simplified sequence diagram of the suggested DA-STE variant. b Schematic extended phase graph of the sequence in (a)

The signal intensities of the pure spin echo (\(_\)) and the isolated stimulated echo (\(_\)) in a double angle configuration with an initial excitation pulse of 90° (as shown in Fig. 1(a)) are given as follows [19, 22,23,24,25]:

$$S_ = M_ \sin \left( _ \cdot \frac} \right)\sin^ \left( _ \cdot \alpha_ } \right)e^ }}$$

(1)

$$S_ = \frac }}\sin \left( _ \cdot \frac} \right)\sin \left( _ \cdot 2\alpha_ } \right)\sin \left( _ \cdot \alpha_ } \right)e^ }} e^ }}$$

(2)

In Equations [1] and [2], the magnitude \(_\) is proportional to the proton density. In case of incomplete magnetization recovery, the steady state longitudinal magnetization before excitation can be used instead [26]. TM and TE are the mixing time and the echo time, respectively (see Fig. 1(a)). T1 is the longitudinal relaxation time and \(_\) refers to the nominal flip angle (as indicated in Fig. 1(a)). Conveniently, B1 refers to the relative \(\left|_^\right|\) which is calculated by the ratio of the actual and nominal flip angle (\(\alpha /_\)). When TM << T1, the cosine of the echo intensity ratio directly yields an estimate of the local flip angle \(\alpha\) [19, 27]:

$$\alpha = \text_ \cdot \alpha_ \approx \cos^ \left( }} }}} \right)$$

(3)

Thus, the sequence provides an estimate of the local flip angle along with the phase from the pure spin echo SE. It is worth noting that Equation [3] is independent of the flip angle of the first excitation pulse as it cancels out in the signal ratio \(_/_\). Setting \(\alpha\) to 90° (as shown in Fig. 1(a)) proves to be beneficial for the signal-to-noise ratio (SNR) for the transceive phase estimation (using both \(_\)) and for the estimation of the local B1 map (as can be verified with Monte Carlo simulations). Note that in case of a non-uniform slice profile, a bias is introduced to Equation [3] that needs to be accounted for to obtain the actual B1. Further details are given in section “Simulations”.

The uncertainty of various B1-mapping methods, such as DA-STE among others, has been investigated in-depth by Pohmann and Scheffler [26]. In the frame of its impact on EPT, the proposed method is in addition compared to the robust but lengthy double angle method using two gradient echo acquisitions [21] (GRE-DAM). The uncertainty was estimated by Monte Carlo simulations and analytically by error propagation according to Pohmann and Scheffler [26]. More detailed information is found in the Supplementary Information.

The uncertainty of the spin echo phase was estimated by the inverse of the SNR of the spin echo magnitude image [28]. Due to the averaging of two phase measurements with opposite read out gradient polarity, the uncertainty of the transceive phase is given by \(_^}=\sqrt\cdot _}^\), where \(SN_\) is the SNR of one of the two magnitude images. \(SN_\) was obtained by the ratio of the average signal magnitude within the object and the standard deviation of the real part of the image in the background (air).

Under the assumption that the local EPs vary slowly, the conductivity \(\sigma\) was reconstructed based on the homogeneous helmholtz (hH) equation [12, 29] yielding:

$$\sigma = \frac \omega }}} \left\ B_^ }}^ }}} \right\} = \frac \omega }}\left( \varphi^ + 2\frac_ \cdot \nabla \varphi^ }}_ }}} \right)$$

(4)

where \(_\) is the magnetic vacuum permeability and \(\omega\) is the Larmor frequency. The conductivity \(\sigma\) is mainly related to the leading term on the right side of Equation [4], which contains only the phase. To avoid additional scan time, the second term containing B1 is thus often neglected yielding the phase-based reconstruction formula [3, 11]:

$$\sigma \approx \frac \varphi^ }} \omega }}$$

(5)

The phase-based approximation can result in errors of 20% at locations of high B1 gradients in the brain at 3 T [11]. Therefore, for an accurate estimation of the tissue conductivity, the B1 information is needed (cf. Equation [4]). Because the transmit phase \(^\) is not directly accessible in MRI, it is approximated by half the measurable transceive phase \(^\) [10, 11]. The latter is the superposition of the transmit and the receive phase. These phases are assumed to be similar for symmetrical objects up to field strengths of 3 T if a quadrature birdcage coil-like configuration is used for transmission and reception. The full hH equation and the magnitude-based expression for the permittivity are given in the Supplementary Information (see Equations [S3] and [S4]).

The local Laplacian and gradients for each voxel were estimated by fitting a second-order polynomial in the voxel’s neighborhood [7]. To mitigate tissue boundary errors, fitting was only applied to voxels that have magnitude values within 25% of the center voxel [7, 12]. The Eps were subsequently smoothed with a tissue boundary preserving median filter [7, 12] using the same strategy. The window size for the Laplace estimation was set to 7 × 7 × 7 voxels and 13 × 13 × 13 voxels for the median filter. Skull stripping and was performed using the standard software package FSL (FMRIB Software Library v6.0, Oxford, United Kingdom) [30]. For segmentation, the TSE brain images were mapped to the standard MNI average brain atlas [31, 32] using FSL’s linear and non-linear registration tool [33,34,35]. The resulting white and gray matter probability masks [36] were registered back to the TSE images, which had the same position and resolution as the DA-STE images. Binary masks were then generated for probability values above 80% to estimate the average conductivity of white and gray matter.

All scans were performed at 3 T (Magnetom Prisma, Siemens Healthineers, Erlangen, Germany) using a dual-tuned 1H/23Na quadrature head coil for transmission and reception (Rapid Biomedical, Rimpar, Germany). Acquisitions were performed both in a phantom and in the brain of a healthy volunteer. In vivo human experiments were approved by the local ethics committee and informed consent was obtained from the volunteer beforehand. The saline phantom contained a salt concentration of 3 g/L, resulting in a theoretically estimated conductivity of 0.5 S/m at the Larmor frequency and at room temperature [37]. Relaxation times were reduced to tissue-comparable relaxation times of T1 ~ 930 ms and T2 ~ 75 ms by adding 0.125 mM MnCl2.

The double angle method using two single gradient echo acquisitions [21] (GRE-DAM) is used as a reference method for the B1 map due to its similarity to the proposed method and its robustness to motion. For GRE-DAM, an estimate of the local flip angle \(\alpha\) is obtained by \(\alpha =}^\left(_/\left(2_\right)\right)\). \(_\) and \(_\) represent the signals of the gradient echo measurements with the flip angles \(_\) and \(_=2_\), respectively. A long TR is used to ensure that the B1 map is largely free from T1 effects. The flip angle of the first and the second acquisition was set to 90° and 180°, respectively, to increase the SNR of the resulting B1 map.

The gold standard of the transceive phase is the average of two spin echo phase acquisitions with opposite readout gradient polarity and was used to verify that the phase from the first echo of the suggested DA-STE sequence indeed corresponds to the conventional spin echo phase and thus to the transceive phase.

The sequence parameters of the two reference methods (SE, GRE-DAM) and the suggested approach (DA-STE) are compiled in Table 1 for the phantom measurements. For DA-STE, the gradient moment A in Fig. 1(a) was set to induce a dephasing of 4π per voxel (~ twice the readout gradient moment) and RF spoiling was enabled for each slice. For the brain measurement, the field-of-view (FOV) was 192 × 168 × 120 mm3.

For tissue-boundary preserving estimation of the Laplacian and median filtering, a T2-weighted image was acquired with a 90°–180° turbo spin echo (TSE) sequence provided by the manufacturer of the system. The TR was set to 10 s with a TE of 105 ms and a receiver bandwidth of 200 Hz/px. The voxel size was 2 × 2 × 2 mm3 and scanning lasted 2:44 min.

All simulations and calculations were performed using MATLAB R2019a (The MathWorks, Natick, MA) unless stated otherwise.

For DA-STE, the \(_\) and \(_\) signals were obtained using a 2D Bloch simulation for an ensemble of 361 × 361 magnetization vectors to simulate the effects within the xz-plane. The effects of the gradient moments in the frequency encoding (x) and slice selection (z) direction are thus accounted for individually. Signal simulation was performed using the parameter settings of Table 1. The impact of selective RF pulses was derived using the hard pulse approximation [38, 39].

Due to selective excitation, the flip angle estimated by Equation [3] must be corrected. To this end, \(_\) and \(_\) signals of DA-STE were simulated for varying relative B1 values between 0.5 and 1.5 with a step size of 0.05. The resulting apparent flip angles (or rather B1) were calculated according to Equation [3] and correction was performed by interpolating the actual (input) B1 as a function of the apparent (output) B1, B1,apparent. For interpolation MATLAB’s interp1 function with the method ‘spline’ was used. Note that the correction was calculated for a fixed T1 and T2 (T1/T2 = 1 s/100 ms) which can result in a potential bias for deviating T1, such as for fluids. The same was done for the GRE-DAM for which T1 effects are smaller due to the longer TR.