Toward accurate and fast velocity quantification with 3D ultrashort TE phase-contrast imaging

Purpose: Traditional phase-contrast MRI is affected by displacement artifacts caused by non-synchronized spatial-and velocity-encoding time points. The resulting inaccurate velocity maps can affect the accuracy of derived hemodynamic parameters. This study proposes and characterizes a 3D radial phase-contrast UTE (PC-UTE) sequence to reduce displacement artifacts. Furthermore, it investigates the displacement of a standard Cartesian flow sequence by utilizing a displacement-free synchronized-single-point-imaging MR sequence (SYNC-SPI) that requires clinically prohibitively long acquisition times. Methods: 3D flow data was acquired at 3T at three different constant flow rates and varying spatial resolutions in a stenotic aorta phantom using the proposed PC-UTE, a Cartesian flow sequence, and a SYNC-SPI sequence as reference. Expected displacement artifacts were calculated from gradient timing waveforms and compared to displacement values measured in the in vitro flow experiments. Results: The PC-UTE sequence reduces displacement and intravoxel dephasing, leading to decreased geometric distortions and signal cancellations in magnitude images,


INTRODUCTION
Phase-contrast (PC) MRI has found widespread clinical use for the quantitative assessment of blood velocities and flow in the diagnosis of various cardiovascular diseases. 1ith the introduction of 4D flow MRI, the velocity vector field in an entire imaging volume can be captured throughout the cardiac cycle. 2,3These velocity measures provide the basis for comprehensive flow analysis and also enable the assessment of additional hemodynamic parameters such as pressure gradients, 4 wall shear stress (WSS), 5 kinetic energy, 6 pulse wave velocity, 7 vortex formation, 8,9 and intracardiac flow component analysis. 10uch parameters can provide additional information for studying disease mechanisms, diagnosis, and treatment monitoring and planning.4D flow MRI is also increasingly used for in vitro imaging of synthetic and biological flow phantoms, for example, to study underlying hemodynamic effects, 11,12 assess accuracy and precision of PC MRI hemodynamic measures, 13 predict surgical outcomes with rapid prototyping models, 14,15 and evaluate computational fluid dynamics modeling results. 16uch in vitro scans offer advantages over in vivo scans, including prolonged scan times to support higher spatial/temporal resolutions and/or higher precision, controlled and repeatable acquisitions without motion artifacts, and comparisons with other velocity measures such as particle image velocimetry. 17,18Whereas several correction methods are routinely applied to PC MRI during image reconstruction and postprocessing 19 to reduce errors due to eddy currents or concomitant fields, some other errors remain unaccounted for in traditional PC MRI.Notably, the non-zero TE causes intravoxel dephasing and associated signal voids; time differences between the spatial-encoding time points of different axes yield geometric distortions or spatial displacements 20 ; and time differences between spatial-and velocity-encoding time points lead to errors in the velocity maps when acceleration is present, [21][22][23][24] sometimes denoted as velocity displacements. 20These artifacts have already been reported in the early development stages of PC MRI 22,23 but have received limited attention since the introduction of 4D flow MRI.However, they become of concern if spatial accuracy is of importance, for example, when (1) comparing MRI velocity fields with other spatially resolved measures such as particle image velocimetry, 17,18 (2) merging data sets such as computational fluid dynamics and PC MRI for improved velocity fields data 25 or segmentations obtained from MR angiography or cardiac anatomy datasets with PC MRI for improved vessel boundary detections, 26 or (3) calculating stream-lines along torturous vessel paths. 27In addition, spatially sensitive measures such as WSS have been shown to be affected. 28ifferent solutions have been proposed to counteract the displacements and geometric distortions in PC MRI.First, displacements can be partially corrected by estimating the actual location of affected voxels by reversing the displacement with a single step 20 or iterative solver 29 retrospectively during data processing.Second, geometric distortions caused by the non-synchronized position-encoding time points 24,30 (t 0 i , for axis i = x, y, z) have been reduced by changing the flow-encoding scheme from the fast but non-synchronous "Bernstein encoding scheme" designed to minimize TE 31 to an echo-encoded scheme, 23,32 in which the displacement is consistent along all three encoding directions, that is, t 0 x = t 0 y = t 0 z .This prospective solution, however, does not reduce the velocity displacement.In contrast, it often rather increases the displacement because of prolonged TEs.Third, the spatial shift as well as the geometric distortions are directly impacted by the choice of imaging parameters that affect the sequence timing parameters, including the chosen velocity-encoding sensitivity (VENC), sequence timing parameters, receiver bandwidth, and spatial resolution.Minimizing TE is not only beneficial for reducing intravoxel dephasing 21,33 but also for reducing the time interval between position-and velocity-encoding ΔT v i along axis i, and thus the velocity displacement. 34nother prospective solution for reducing geometric distortions, spatial and velocity displacements, is to change the type of sequence trajectory.Bruschewski et al. presented a SYNC-SPI (synchronized-single-pointimaging) sequence 35 that eliminates geometric distortions and displacements by simultaneously encoding the velocity and position on all axes.However, whereas the SYNC-SPI sequence provides highly accurate velocity vector fields, it acquires a single k-space point per repetition time, yielding extraordinarily long scan times of several hours up to several days that are often highly demanding concerning the stability of the flow setup and MRI system or even unfeasible.In practice, a faster method than the SYNC-SPI is often needed, but at the same time it should provide more accurate results than conventional Cartesian acquisitions and independence of the chosen orientation.
Using asymmetric Cartesian readout sampling schemes 36 can reduce ΔT v i , but still the artifacts are nonisotropic.Further shortening of ΔT v i can be achieved by readouts that start in the center of acquisition k-space, for example stack-of-stars 37,38 or stack-of-spirals. 34,39It was shown that these sequences are beneficial for investigating complex flow in stenotic jets with high velocities 34,38,39 since short TE values reduce intravoxel dephasing and signal loss.
For magnitude-based applications, UTE sequences 40 combined with nonselective RF pulses have been proposed that use isotropic center-out acquisitions, for example, "kooshball" 40 -like radial trajectories.Such radial UTE trajectories have also been used in combination with bipolar gradients for PC imaging, which has been termed PC-UTE.We adopt the term PC-UTE in this work, although an ultrashort TE value cannot be realized due to the flow-encoding gradient duration prior to the readout.In other works, 3D radial PC-UTE has been investigated in animal studies: improved velocity mapping compared to a corresponding Cartesian sequence was shown for self-gated 4D PC-UTE 41 with a delayed readout gradient applied after the flow-encoding gradients.WSS values were derived from a self-gated 4D radial PC sequence 42 utilizing an asymmetric readout.Further, peak flow and maximal velocities were obtained from accelerated 4D radial "center-out" PC data 43 utilizing compressed sensing with both readout gradients and data sampling delayed.Most of these works focused on the beneficial performance of radial center-out PC sequences regarding flow artifacts, 41,43 robustness to motion, 43 and intravoxel dephasing. 41However, a thorough investigation of the associated displacement artifacts is lacking.
In this work, we investigate a 3D PC-UTE sequence and quantify its displacement artifact in an in vitro flow experiment.Comparison with a Cartesian PC sequence is performed, and SYNC-SPI data is acquired to serve as the ground truth.ΔT v i -values are calculated from the sequence timings for different spatial resolutions and image orientations, and displacement predictions from these calculated ΔT v i -values are compared to the measured values.Results obtained by the 3D PC-UTE sequence are compared to a Cartesian acquisition post-processed with the method proposed by Thunberg et al. 29 to correct the displacement artifacts.

Calculation of encoding time points and displacement values
Although the position-and velocity-encoding process in PC MRI occurs within a finite duration, both processes are often considered to happen at single time points. 24,30,31or accurate flow encoding, the placement of these time points, that is, (t 0 x , t 0 y , t 0 z ) for encoding the position along the x-, y-, and z-axis, as well as (t 1 x , t 1 y , t 1 z ) for encoding the velocity along the axes, are essential. 28,35,31,24They depend on various factors, including the sequence trajectory, image orientation, different (velocity) encoding schemes, gradient performance, and various sequence parameters such as receiver bandwidth, spatial resolution, and TE.
Along the readout (RO) direction 30 of a Cartesian sequence, for example, t 0 x is considered to be located at the center of the echo and coincides with the center of the readout gradient for a symmetric readout (Figure 1).In contrast, if the fastest velocity encoding is implemented according to Bernstein et al., 31 t 0 y , t 0 z along the phase encoding (PE) directions are located at the "gravity time point" of the applied PE gradients. 24Thus, the time points t 0 i are not synchronized between the axes, not only in the Cartesian case but also for other acquisitions such as stack-of-stars.As a result of this non-zero time interval ΔT 0 ij = t 0 i − t 0  between position-encoding 24 and axes ij = xy, xz, yz, the position of moving spins appears spatially dislocated and can lead to geometric distortions.In PC MRI, the time points for encoding the velocity are at the "gravity time point" of the bipolar gradient difference function. 24In most sequences, the bipolar gradients are switched simultaneously (although potentially with different zeroth moments m 0 and first moments m 1 ) along the three axes; thus, the time points t 1 x , t 1 y , t 1 z are typically close in time (Figure 1).However, there are time differences ΔT v i between t 1 i and t 0 i that can vary between the axes, with the readout axis typically showing the largest ΔT v RO .Since the velocity is encoded prior to the spatial position, the peak velocity in a stenotic region is shifted downstream.This velocity displacement Δd i of moving spins with velocity component v i and time interval ΔT v i along axis i is given by To investigate the dependency of the displacement artifacts on the spatial resolution and associated encoding time points, the Cartesian sequence is simulated with the vendor's integrated development environment for applications (Siemens Healthineers, Erlangen, Germany, IDEA, version VB17) for isotropic resolutions ranging between 0.6 and 5.3 mm (TE = 5.4-2.8ms), with otherwise identical settings (VENC: 2.5 m/s, receiver bandwidth: 450 Hz/Px, maximal gradient amplitude: 19.2 mT/m, and maximal gradient slew rate: 96 mT/m/ms).In the PC-UTE sequence, the time points t 1 i and t 0 i , as well as ΔT 0 ij and ΔT v i , are not affected by changing the resolution or receiver bandwidth because the "center-out" sampling requires no m 0 for prewinding (Figure 1).Therefore, the radial sequence is simulated only once for an applied VENC of 2.5 m/s.From the obtained gradient waveforms, the encoding timepoints (t 1 i , t 0 i ) and ΔT v i are derived and the expected velocity displacement values Δd i calculated (Table 1).

F I G U R E 1
Sequence diagrams for flow acquisitions, exemplarily for one TR.From left to right: SYNC-SPI (reference), 35 Cartesian, 31 and proposed 3D PC-UTE sequence.In Cart RO ∥ , RO, PE 1 , PE 2 are played out on z-, x-and y-axis (RO parallel to z-axis).In Cart RO ⊥ , RO, PE 1 , PE 2 are played out on x-, z-, and y-axis (RO orthogonal to z-axis).Position-encoding time points (t 0 x , t 0 y , t 0 z ) are marked by black dashed lines, and velocity-encoding time points (t 1 x , t 1 y , t 1 z ) by gray dotted lines.Bottom row: a sketch of acquired data.PC, phase contrast; PE, phase encoding; RO, readout; SYNC-SPI, synchronized-single-point-imaging.

Phase-contrast MRI sequences
In this work, three sequences as outlined in the following are used to investigate the displacement artifacts quantitatively.

3D Cartesian PC sequence
A Cartesian 3D PC MRI sequence is investigated, including different imaging orientations.The velocity-encoding gradients are implemented according to the established Calculated and measured displacement values Δd z with corresponding time intervals ΔT v z of encoding directions played out along stenosis for flow experiments with constant flow rates Q 1 ∕Q 2 ∕Q 3 = 40∕30∕20 mL/s (corresponding to maximal expected velocities v z = 2.0/1.5/1.0 m/s), VENC = 2.5 m/s, bandwidth = 450 Hz/Px, 1 mm isotropic resolution.method proposed by Bernstein et al., 31 providing the fastest encoding scheme and a minimal TE. 31 Here, the m 0 of the pre-/rewinding gradients are included in the bipolar gradients for velocity encoding, yielding non-synchronized position-and velocity-encoding time points for the three individual axes (Figure 1).Consequently, ΔT 0 ij and ΔT v i vary between axes resulting in both velocity displacements and geometrica distortions.

3D PC-UTE sequence
3D velocity-encoding gradients with the shortest possible duration 31 are integrated into a 3D UTE radial readout sequence 40 in which the readout starts in k-space center and subsequently moves radially outwards.A rectangular-shaped 500 μs-long RF pulse is followed by bipolar gradients (total duration = 940 μs for VENC = 2.5 m/s, TE = 1.2 ms), trapezoidal readout, and spoiling gradients (Figure 1).The first readout-spoke ends at the south pole (z-axis) and the ends of subsequent spokes spiral up to the north pole, while each TR is repeated four times with different velocity encodings.While the readout gradient changes direction, the four velocity-encoding gradients remain constant along the axes throughout the sequence.Due to "center-out" acquisitions, the m 0 of the bipolar-encoding gradients vanishes.A two-sided velocity-encoding scheme 31 is implemented with a constant absolute m 1 across all four flow encodes (FE) and all three directions (x, y, z) but with swapping signs: ) , ) .As a result, the velocity-encoding timepoints t 1 i are synchronized for all three axes i = x, y, z.Additionally, the position-encoding timepoints t 0 i are likewise synchronized, which eliminates geometric distortions, and are located at the beginning of the readout.Thus, this sequence shows a minimal time interval ΔT v i along all axes, resulting in an isotropic sensitivity to spatial displacements.

SYNC-SPI
A SYNC-SPI sequence 35 is applied in addition to the other two acquisitions, which is regarded as a ground truth here.The SYNC-SPI achieves ΔT 0 xy,xz,yx = ΔT v x,y,z = 0 for all three axes by synchronizing t 1 i and t 0 i , resulting in a single time point, that is, t 0 x = t 0 y = t 0 z = t 1 x = t 1 y = t 1 z , preventing velocity displacements and geometric distortions.This is feasible because only a single point in k-space is acquired per TR, and the k-space position does not change Sketch of stenotic flow phantom setup embedded in polyvinylpyrrolidone and positioned along physical scanner coordinates (x, y, z).The blue arrow marks the inflow direction.
during readout but it comes at the expense of approximately two orders of magnitude longer acquisition time, that is, multiple hours instead of a few minutes.

Setup and flow measurements
In vitro flow experiments are performed on a 3T MR system (Siemens Verio, Erlangen, Germany) with a flow pump (CardioFlow 5000 MR, Shelley Medical Imaging Technologies, London, Canada) using three different constant flow rates Q 1 /Q 2 /Q 3 of 40/30/20 mL/s leading to maximal expected velocities of 2.0/1.5/1.0 m/s within the stenotic region of the flow phantom.Distilled water is pumped through the phantom consisting of a Venturi nozzle followed by a U-bend of 150 mm diameter (Figure 2) embedded in static 1.5% agarose gel, representing a stenotic aorta. 28The flow phantom has an aortic diameter of 9.5 and 15.0 mm proximal and distal to the stenosis (Figure 2).Reynold numbers of approximately 2700/2000/1300 and 5100/3800/2500 are expected in the pipe inflow and the stenosis for All 3D flow sequences are acquired with nontime-resolved three-directional velocity encoding (3D PC MRI) using a 15-channel knee coil (QED, Mayfield Village, OH USA).For characterization of the MR gradient system and subsequent correction of the data due to deviations from the nominal radial readout trajectory, the gradient impulse response function (GIRF) 44 is measured using 12 triangular gradients (durations of 50-160 μs) with a magnetic field camera (Skope Magnetic Resonance Technologies, Zurich, Switzerland) according to Vannesjo et al. 44 The flow phantom experiments are performed with the Cartesian, PC-UTE, and SYNC-SPI sequences for Q 1 , Q 2 , and Q 3 .While ΔT 0 ij and ΔT v i are independent of the image orientation for the PC-UTE and SYNC-SPI, this is not the case for the Cartesian sequence.Hence, the Cartesian acquisitions were performed twice, first with readout-direction parallel (Cart RO ∥ ) and second with readout-direction perpendicular (Cart RO ⊥ ) to the stenosis.Figures 1 and 2 display all sequences for a single TR, with corresponding encoding time points and the phantom setup/geometry.
All flow sequences are acquired with a nominal 1 mm isotropic resolution and 12

Reconstruction and post-processing of velocity data
All datasets are reconstructed using custom-written reconstruction software in Python (version 3.7).To correct for gradient imperfections in the radial acquisitions, a GIRF 45 is applied to both the nominal bipolar and radial readout gradients.The radial datasets are then reconstructed with density compensation, and a nonuniform fast Fourier operator is applied to the GIRF-predicted k-space trajectory.
No Maxwell correction is applied throughout this work.Background phase correction using a third-order polynomial fit on static material (see Figure 2 for surrounding tissue) is applied for each velocity dataset, and the resulting fit is subtracted from the velocity maps.

2.5
Data evaluation

Velocity displacements
The highest velocity shifts within the phantom are expected in the stenotic region.To quantify the shift and its dependence on the image orientation, the velocity component v z along the stenosis is evaluated for Q 1 − Q 3 and all three acquisitions.For this, the v z -component of the SYNC-SPI data is shifted by v z ⋅ ΔT v z for different ΔT v z -values.Since the expected displacement values are smaller than the acquired voxel size, the velocity data is first interpolated to a finer grid, shifted by amount ΔT v z , and then interpolated to the original voxel grid again.Subsequently, the correlation between the shifted SYNC-SPI and the other acquisitions is obtained.The ΔT v z with the best correlation to the measured v z of Cart RO ∥ , Cart RO ⊥ , and PC-UTE is determined, and the measured velocity shift Δd z for the velocities of 2.0/1.5/1.0 m/s is derived.

Geometric distortions
Geometric distortions are present in both the magnitude and velocity maps due to non-zero and non-synchronized ΔT 0 ij -values, which may result in the signal of moving spins being shifted into the walls of the vessel pipes.To visualize the geometric distortions in the velocity data, masks of the vessel boundaries are overlaid with the velocity maps for the different PC MRI acquisitions.Masks are obtained with the software ITK-SNAP (version 4.0.0)by applying seed points in the magnitude images for each acquisition.The shift of the fluid within the pipe is determined by comparing the masks to the SYNC-SPI-mask and the maximum shift distance is reported.

Sensitivity to intravoxel dephasing
For acquisitions with high TE, intravoxel dephasing can lead to signal cancellations and dropouts, which may affect velocity quantification and derived parameters.To investigate the performance of the Cartesian, PC-UTE, and SYNC-SPI sequences regarding signal attenuations, signal intensities in the magnitude data are analyzed in a region of interest (ROI) of 16 voxels downstream of the stenosis, which present the largest decrease of signal intensities due to complex flow for Q 1 − Q 3 .For all acquisitions, the signal intensities are normalized to the mean inflow signal prior to the stenosis.The largest signal dropouts downstream of the stenosis are compared to those of the SYNC-SPI data.

Flow and velocity value evaluation
To investigate the performance of all flow sequences, both flow rate and velocity component values are investigated.
First, the mean velocities with SDs are determined within ROIs of 36/122/11 voxels at different cross-sectional areas (pipe inflow/outflow/stenosis) of the flow phantom and Q 1 .Volumetric flow rates are then derived by multiplying these areas with the mean velocities.The mean velocity component v x with SD is determined in a ROI of 43 voxels for all acquisitions in the arc of the pipe.Since the radial data is more strongly affected by system imperfections compared to the Cartesian sequence, the velocity values of v z along a velocity line are investigated for the PC-UTE sequence regarding the performance of GIRF and background phase corrections.

Retrospective displacement correction method
A retrospective correction method proposed by Thunberg et al. 29 is implemented and applied to both the velocity and magnitude data obtained with Cart RO ∥ and flow rate Q 1 , for which the displacement is most prominent along the z-axis.As described by Thunberg et al., the velocity values in the acquired data are set back by the displaced distance given by the velocity values multiplied by the time interval ΔT v RO .From this corrected velocity field, streamlines are calculated and tracked over ΔT v RO .The displaced velocities are then shifted backward from the end to the starting position of the streamlines, resulting in a further corrected velocity field.This process is repeated iteratively several times, including streamline calculation based on the new corrected velocity fields.

Calculation of encoding time points and displacement values
In Figure 3A-C, bipolar gradients with corresponding encoding time points and time intervals (shaded region) are displayed for one flow encode for the Cartesian and PC-UTE sequence as used in the measurements for a 1 mm isotropic resolution.The Cartesian encoding resulted in predicted values ΔT v z of 3.3 ms for the readout and 0.7 ms for the first PE direction.For the PC-UTE sequence, an isotropic value of ΔT v i = 0.5 ms was calculated for all three spatial axes.With these time intervals, displacement values Δd z of 6.6/5.0/3.3 mm for Cart RO ∥ , 1.4/1.0/0.7 mm for Cart RO ⊥ , and 1.0/0.8/0.5 mm for PC-UTE were calculated for the three maximum velocities of Q 1 ∕Q 2 ∕Q 3 , respectively.Calculated displacement values and corresponding time intervals are listed in Table 1.
Figure 3D-F illustrates the absolute displacement in mm as a function of the resolution, which was calculated for the three different sequences with VENC = 2.5 m/s.The SYNC-SPI shows no displacement independently of the resolution and for all axes.For the Cartesian sequence, the displacement depends on the spatial resolution, encoding axis, and velocity direction.The displacement along Cartesian readout dominates for a velocity of 2.0 m/s, increasing from 5.5 to 8.1 mm for isotropic resolutions of 2.0 and 0.6 mm, respectively.This effect becomes even more apparent when calculating the shift relative to the voxel size, as shown in Figure 3G-I.Here, the displacement increases from 2.7 to 13.4 voxels.For the PC-UTE, the calculated displacement is independent of both the direction and resolution but shows a non-zero value of 1.0 mm for a velocity of 2.0 m/s.

Flow measurements
For the four flow acquisitions, the magnitude maps for Q 1 (Figure 4A) and velocity maps for the z-component of all three flowrates Q 1 − Q 3 (Figure 4B-D) are displayed in Figure 4.For each flow rate, all three velocity components and magnitude maps are displayed in Figures S1-S4 for the same slices.The displacement shifts and geometric distortions are less prominent with smaller flow rates and correspondingly smaller velocity values.

Velocity displacements
The largest displacement occurs in Cart RO ∥ since the readout axis is aligned along the highest velocity component v z (marked by gray arrow in Figure 4).For Q 1 − Q 3 , the velocities along the dashed white line (Figure 5B) are displayed in Figure 6 to investigate the velocity shifts.Measured velocity displacement Δd z resulted in 6.0/5.0/3.0 mm for Cart RO ∥ , in 1.0/1.0/1.0 mm for Cart RO ⊥ , and in 1.0/1.0/1.0 mm for PC-UTE compared to the SYNC-SPI acquisition for Q 1 /Q 2 ∕Q 3 , respectively.Cart RO ⊥ and PC-UTE show similar results for Q 1 − Q 3 (Figure 6).In Figure S5, slices with sagittal orientations show similar results compared to slices with coronal orientations (Figure 4), depicting the dependence of the artifact in 3D volume.The dependence of the displacement on the acquired resolution of Cart RO ∥ and Cart RO ⊥ is illustrated in Figures S6 and S7.Velocity shifts of Cart RO ∥ compared to Cart RO ⊥ resulted in 9/2/2/1/1 voxels corresponding to 5.4/3.2/3.6/2.3/2.7 mm for the isotropic resolutions 0.6/1.6/1.8/2.3/2.7 mm.

F I G U R E 3
Gradient waveforms of bipolar gradients for Cartesian RO direction (A), Cartesian PE 1 -axis (B), and PC-UTE sequence (C) (VENC = 2.5 m/s, bandwidth = 450 Hz/Px, 1 mm isotropic resolution).Readout gradients are displayed in (A) and (C).The position-encoding time points (t 0 z ) are marked by black dashed lines, the velocity-encoding time points (t 1 z ) by gray dotted lines.The gray shaded regions mark the time interval ΔT v z between t 1 z and t 0 z .Dependence of velocity displacement on isotropic resolution (range: 0.6-5.3mm) in mm (D-F) and in number of voxels (G-I) for a flow velocity of 2.0/1.5/1.0 m/s.The larger displacement for Cart PE 1 for higher resolutions is caused by increasing area and length of the bipolar gradients along PE-axis with decreasing resolution; t 0 z (located at "gravity time point" 24 of the second lobe) is more affected than t 1 z (located at "gravity time point" of bipolar gradient difference, that is, more toward the beginning of second lobe).All flow sequences were calculated with the vendor's integrated development environment for applications.VENC, velocity-encoding sensitivity.

Geometric distortions
In Figure 4, the SYNC-SPI and PC-UTE data show no geometric distortions.However, in both Cartesian acquisitions, geometric distortions are observed.White arrows pointing to the pipe bending indicate the dependency of the distortions on the imaging axes in Figure 4.
For Cart RO ∥ , with the readout direction applied along the bottom-top direction, the signal of the moving fluid within the pipe is shifted 2-3 mm upward such that the wall appears to be thinner.Since the velocity component v z shows a range of different values, moving spins are shifted by varying amounts.If the readout direction is oriented along the left-right direction (Cart RO ⊥ ), the moving fluid is shifted toward the right within the arch, and the pipe wall appears 2-3 mm broader compared to the SYNC-SPI.
F I G U R E 4 3D flow measurements in a stenotic flow phantom with VENC = 2.5 m/s and flow rates Q 1 − Q 3 .(A) displays magnitude images of Q 1 , (B-D) the velocity maps of the velocity along the stenosis (v z ) for each flow rate.White arrows mark geometric distortions in both magnitude and velocity maps.Gray arrow mark velocity displacement mostly present in the velocity map v z of Cart RO ∥ .Slices 9, 20, 38, 120 are displayed for SYNC-SPI, Cart RO ∥ , Cart RO ⊥ , and PC-UTE, respectively.

Sensitivity to intravoxel dephasing
Intravoxel dephasing can lead to signal cancellations in both velocity and magnitude maps.Whereas the SYNC-SPI shows the fewest signal fluctuations, both Cartesian acquisitions show a substantial signal loss in the magnitude images (black regions in Figure 5A), which correspond to complex flow regions and can be reduced by shortening TE.Signal intensities evaluated along the dashed white line in Figure 5A are displayed in Figure 7 for The signal intensity decreased maximally by 92%/64%/23% (Q 1 ), 78%/42%/19% (Q 2 ), and 51%/25%/20% (Q 3 ) for Cart RO ∥ /Cart RO ⊥ /PC-UTE compared to the SYNC-SPI.Signal dropouts do not only become apparent in the coronal view shown in Figure S4 but also in magnitude maps with sagittal orientations as depicted in Figure S5.
Additionally, the dependence of signal cancellations for different resolutions, and thus TE values, is illustrated in Figure S6.For the highest resolution (0.6 mm 3 ), signal cancellations are mostly prominent and decrease for lower resolutions (>1.0 mm 3 ).

Flow and velocity value evaluation
Mean velocity and volumetric flow rates are obtained for different cross-sections of the flow phantom for Q 1 (Figure S8).In the region of the stenosis, volumetric flow rates are underestimated by 47%/5%/10% for Cart RO ∥ ∕Cart RO ⊥ ∕PC-UTE compared to the SYNC-SPI.The high underestimation in Cart RO ∥ demonstrates the velocity displacement of high velocities shifted downstream of the stenosis.The velocity component v x

F I G U R E 5
Zoomed magnitude (A) and velocity maps (B) downstream of the stenosis exemplarily for flow rate Q 1 .Signal cancellations in (A) are evaluated along the white dashed line in Figure 7. Velocity displacement shifts of highest velocities marked by white lines in (B), evaluation of velocities along the white dashed line in Figure 6.The velocity colormap corresponds to Figure 4.

F I G U R E 6
Quantification of velocity displacements in (A-C) corresponding to flow rates Q 1 − Q 3 .Velocity components v z are plotted along the dashed white line in Figure 5

Retrospective displacement correction method
The velocity and magnitude maps obtained after applying the correction method 29  are displayed in Figure 8. Displaced velocities (v z ) are successfully shifted backward (Figure 8C).The resulting shift with Cart RO ∥ obtained after correction deviates less than 1 mm to Cart RO ⊥ , PC-UTE, as well as SYNC-SPI.However, the corrected velocity values do not completely match with the SYNC-SPI data.The signal of displaced spins is successfully shifted backward.Nevertheless, this correction method does not remove signal cancellations due to intravoxel dephasing.

DISCUSSION
In this work, a 3D PC-UTE flow sequence was developed and investigated in in vitro flow experiments with a focus on displacement artifacts and intravoxel dephasing.To the best of our knowledge, it is the first report in the literature that quantitatively analyzes the displacement artifact with a PC-UTE sequence and compares it to both the displacement of a Cartesian and ground truth SYNC-SPI sequence.
Over the past decade, the increasing availability of 4D flow sequences, post-processing packages, and 3D printers has led to a rising number of in vitro PC flow imaging studies investigating flow patterns and hemodynamic parameters in replicas of healthy and diseased vessels in detail.This approach enables analyzing the impact of specific anatomic features on the underlying flow pattern.It enables prolonged scan times not feasible in in vivo studies, thus supporting higher spatial resolutions.However, applying Cartesian-based, PC flow sequences or PC sequences with long TE values for high-resolution applications yields strong spatial and/or velocity displacement artifacts, 20,28 which can vary between the readout and PE axes. 22,28These effects impact spatial accuracy as well as measured velocity amplitudes, velocity vector fields, and derived hemodynamic parameters such as the WSS 28 and kinetic energy.
For the proposed PC-UTE sequence, ΔT v i only depends on the VENC value and gradient performance.Thus, newer MR systems with higher gradient performance are favorable for achieving shorter ΔT v i durations.In addition, the resolution and the receiver bandwidth impact the displacement for Cartesian sequences.For many flow applications, particularly for WSS estimation, high resolutions are required for obtaining accurate values. 46Relative to the voxel size, the displacement increases disproportionally high, as shown in Figure 3E-I, demonstrating the importance of sequences that minimize those effects when high spatial resolutions are indispensable.As shown by Schmidt et al., 28 varying WSS values are obtained from velocity fields acquired with either "Bernstein encoding scheme" 28,31 (t 0 x ≠ t 0 y , t 0 z ) or a scheme that synchronizes t 0 i for all axes at TE 28 (t 0 x = t 0 y = t 0 z ).Since the velocity-encoding scheme in our PC-UTE sequence eliminates ΔT 0 ij and minimizes ΔT v i , WSS values derived from this acquisition are expected to be closer to the ground truth due to smaller and isotropic velocity shifts.
Synchronized encoding and vanishing ΔT 0 ij and ΔT v i values are achievable by the SPRITE 47 or SYNC-SPI   method, but they come at the cost of extensively long scan times.Different types of multi-slice 2D PC-UTE 34,38,39 acquisitions have been investigated for flow imaging and can serve as a tradeoff between reduced scan time over point-based methods and greater accuracy over Cartesian acquisitions.
The PC-UTE scan time is substantially shorter than the SYNC-SPI counterpart but still requires around 25 min in the present implementation, which is 38% longer than the Cartesian scan time.Extending the technique toward 4D flow, including temporal dimensions, is straightforward but further prolongs the scan time.Nevertheless, radial acquisition schemes are well suitable for subsampling 48 and iterative reconstructions, 43 which have not been applied here.
Importantly, UTE sequences also reduce intravoxel dephasing, which was particularly beneficial in complex flow. 39,49In this regard, center-out sampling strategies such as a center-out radial "kooshball" 41,40,50 or cones 51 are preferred for 3D PC imaging.Magnitude and velocity maps show fewer signal cancellations in stenotic regions due to shorter TE, which can be beneficial for quantification of turbulences.
Non-Cartesian UTE acquisitions, however, are susceptible to system imperfections. 52Gradient delays and trajectory deviations 52 must be carefully considered and corrected during reconstruction.We chose the GIRF prediction model 44 to correct for gradient deviations.Correcting for both bipolar and readout gradients included the different polarities of the flow-encoding gradients in front of the radial readout for each flow encode.Eddy currents from the flow-encoding gradients that influence the readout could be mitigated, for example, by delaying the readout; but this would come at the expense of a prolonged ΔT v i .Further investigations need to be carried out on the performance of the applied GIRF model, for example, by measuring the gradients directly with a field camera.
Measurements with varying flow rates demonstrated the linear dependency of the displacement on the velocity values.Consequently, a lower flow rate results only in moderate displacements.Thus, for pulsatile flow, which was not investigated in this work, only minor differences in the velocity displacements between the proposed PC-UTE and Cart RO ∥ are expected for the diastolic phase, whereas larger differences are expected in the systolic phase.Moreover, applying a smaller VENC value for better velocity-to-noise ratio typically increases TE and the length of the bipolar velocity-encoding gradients, thus the risk for geometric distortions and velocity displacements.Moreover, the artifact depends on the underlying geometry; for example, velocity shifts are visible in regions with acceleration and the shift occurs in the 3D volume.Note that geometric distortions also exist in non-PC MRI sequences, which can become important if images from different sequences are merged, for example, if a vessel segmentation for PC MRI is performed on a sequence different from the PC MRI sequence.
Applying the correction method 29 on the Cartesian data successfully shifted high velocities backward.However, signal dropouts that occur in regions downstream of the stenosis could not be retrieved.For the future, further improvements might be achieved by applying a mask of static tissue, preventing the signal from being shifted into the velocity fields.

CONCLUSION
Due to the shortened and directionally independent time interval between spatial-and velocity-encoding, PC-UTE sequences are highly suitable for measuring the underlying velocity vector field more accurately compared to the Cartesian counterpart and in shorter scan times compared to sequences without displacement artifacts.This benefit becomes especially important when high resolution is desired or when imaging high flow velocities.The proposed 3D PC-UTE sequence may be used for detailed investigations of biomedical flows, for example, when studying flow patterns in a 3D-printed replica of vessels, aneurysms, or artificial heart valves.
for four acquisitions.Black arrows mark the large displacement of Cart RO ∥ acquisition (RO; Cart RO ⊥ and PC-UTE show similar performance.was evaluated for all four acquisitions in the bending of the pipe and resulted in −0.20 ∓ 0.01/−0.26∓ 0.01/ −0.23 ∓ 0.01/−0.22∓ 0.01 m/s (Q 1 ), in −0.21 ∓ 0.01/ −0.20 ∓ 0.01/−0.15∓ 0.01/−0.17∓ 0.01 m/s (Q 2 ), and in −0.13 ∓ 0.01/−0.11∓ 0.01/−0.12∓ 0.01/−0.14∓ 0.01 m/s (Q 3 ) for SYNC-SPI/Cart RO ∥ /Cart RO ⊥ /PC-UTE, respectively.The velocity component v y of PC-UTE shows an overall similar velocity pattern to the SYNC-SPI data (Figures S1-S3), whereas Cart RO ∥ and Cart RO ⊥ show deviations in the pipe-bending region, which may result from nonsynchronized position-encoding time points in 3D volume.The velocity component v z of the PC-UTE acquisition is displayed for the nominal and the GIRF-corrected trajectories as well as after background phase correction in Figure S9.
on data acquired with Cart RO ∥ F I G U R E 7 Quantification of signal intensities from magnitude maps for SYNC-SPI, Cart RO ∥ , Cart RO ⊥ , and PC-UTE in (A-C) corresponding to flow rates Q 1 − Q 3 .Signal intensities are plotted along the dashed white line in Figure 5; black arrows mark signal dropouts; and dotted black lines mark pipe edges.

F I G U R E 8
Retrospective correction 29 of velocity data for Cart RO ∥ acquisition ("Cart RO ∥ : corr") applied on (A) velocity component along the stenosis (v z ) and (B) magnitude for flow rate Q 1 , VENC = 2.5 m/s.In (C), velocities along the dashed line in (A) are displayed.Slices 9 and 20 are displayed for SYNC-SPI and Cart RO ∥ .

Table S1 .
Acquisition parameters of 3D flow measurements with constant flow rates Q 1 − Q 3 .Number of slices correspond to PE 2 -steps in Cartesian and SYNC-SPI and to base resolution in PC-UTE acquisition: 26, 60, 60, 256 slices for SYNC-SPI, Cart RO ∥ , Cart RO ⊥ and PC-UTE, respectively.