A Statistical Robust Approach to Design Parallel Transmit Radiofrequency Excitations in MRI

,


Introduction
In magnetic resonance imaging at high magnetic field, the apparition of standing wave effects [1][2][3] prevents uniform excitation of the spins across extended anatomical regions. For brain imaging, these effects are clearly visible at 3 tesla (T) and are often deemed unacceptable at field strength > 4 T. Parallel transmission (pTX) is a promising technology to mitigate this problem as it provides degrees of freedom to "shim" the RF field distribution [4] (RF shimming). e transmit SENSE technique [5][6][7], enabled by the pTX technology, furthermore allows reducing the duration of multidimensional pulses [8] while preserving spatial definition. It thus offers one more solution to tackle the transmit field heterogeneity problem. e rationale here is to regard the problem of achieving a uniform magnetization profile as one particular problem that multidimensional RF pulses can solve, whereby the target magnetization profile is spatially uniform. e design of such pulses, here referred to as flip angle (FA) shimming, plays now an important role at UHF magnetic resonance to enable uniform spin excitations.
us, the so-called fast-kz or spokes pulses [9,10] and the k T -point pulses [11] are well-known slice-selective and nonselective FA shimming strategies taking advantage of the parallel transmission technology to produce highly uniform excitation profiles at a reasonable cost in terms of pulse duration as compared to RF shimming.
However, in addition to being nonuniform in space, the transmit RF field (B + 1 ) also depends on the geometry and electric properties (dielectric permittivity and conductivity) of the sample and thus can vary significantly across individuals. us far, the usual technique to cope with this is to characterize the actual B + 1 profiles using MRI-based measurements and to tailor RF shimming to this particular TX field distribution. is obviously can be time-consuming, in particular when the number of TX channels is high [12]. One alternative approach to handle the variability of the TX field distribution, called universal pulse (UP) [13], is to replace the subject-specific TX field distribution with a collection of such TX field distributions acquired on a database of subjects. e rationale here is to achieve, simultaneously across all subjects of the database, good excitation performances, and assume that this property generalizes to the population that the database is supposed to cover (e.g., the adult population). e validity of this approach has been demonstrated experimentally for various applications and hardware setups [14][15][16]. In those studies, a database of 10 to 20 individuals was deemed sufficient to yield good performances. However, the amount of data remained too limited to analyze quantitatively how pulse performance changes with the design database, and in particular with its size. Nevertheless, it is reasonable to assume that UP performance increases with the size of the design database; conversely, if a too small database is used to compute UPs, there is a risk of "data overfitting," i.e., that the solution displays a higher FA shimming performance on the design database than its actual performance on the entire population.
To better control the actual performance of UP on the entire population, we propose in this work to model the transmit RF field as a stochastic process, whereby every new subject (giving rise to a new B + 1 map) is a realization of this process. is representation appears very useful as it allows formulating the action of a FA shimming pulse in statistical terms. From the RF pulse design perspective, it also provides an adequate framework to apply the socalled statistical robust design theory, proposed by Alotto et al. [17], applied to the design of electromagnetic devices. For the latter application indeed, manufacturing errors of the components or uncontrollable environmental conditions such as temperature or humidity can lead to uncertainties in the response of the system. Applying classical optimization to such a problem can severely degrade the actual performance of the solution. A robust design approach mitigates this problem by accounting for the randomness of the response of the system in the optimization process using an appropriate statistical model. is idea can naturally be extended to the RF pulse engineering problem.
We thus present in this work a statistical representation of the TX field distribution B + 1 in which, for every position r, B + 1 (r) is assimilated to a Gaussian random process. We then derive a method to apply the aforementioned statistical robust design method to the design of FA shimming UPs, which uses the second-order statistics of B + 1 (r) instead of the collection of B + 1 maps as input parameter. An experimental validation of the method is provided with data acquired with an eightfold transmit RF coil for 7 tesla head MRI. We demonstrate in particular (1) the utility of statistical robust design approach to boost computational efficiency, (2) the capacity of the method to enable the integration of a large database TX field maps, and (3) the utility of a statistical framework to compute flip angle statistics and to derive simple performance metrics for quality assurance.

Statistical Modelling of the Resonance TX Field. Let B +
1,c (r) (tesla/volt) denote the spatially dependent (r referring to the position) resonant transmit RF magnetic field of the c th resonator of the N c -fold transmit array and c � 42.57 × 10 6 Hz/T the gyromagnetic ratio of the proton. e so-called control field components ω 1,j � 2πcB + 1,j (rad/ s), 1 ≤ j ≤ N c , may then be arranged in a vector ω 1 (r) � (ω 1,1 , . . . , ω 1,N c ) t ∈ C N c ×1 called the control field vector. e control field vector being subject dependent is viewed as a stochastic variable where each new subject gives rise to a realization of ω 1 . For every location, we assume that the mean value ω 1 (r) and the covariance matrix ( * for the Hermitian conjugate) (ω 1 (r) − ω 1 (r) )(ω 1 (r)− ω 1 (r) ) * 〉 of the control field vector is well defined and we denote them, respectively, by μ(r) ∈ C N c ×1 and C(r) ∈ C N c ×N c . In the following, the control field ω 1 (r) is assumed to have Gaussian statistics, i.e., ω 1 (r) ∼ N(μ(r), C(r)). Let be N independent control field distributions (in practice TX field maps measured in N different subjects). If N is sufficiently large, this dataset can be exploited to estimate the mean μ(r) and covariance matrix C(r) of the underlying random process by computing the sample mean and covariance matrix for this ensemble: For the covariance matrix, we shall use the Ledoit-Wolf [18] correction of the sample covariance estimator when the condition N ≫ N 2 c (N 2 c � 64 in this study) is not satisfied. Given a statistical model for B + 1 , we now show that this information is sufficient to compute universal FA shimming pulses and to study the robustness of such pulses with respect to the intersubject variability of B + 1 . We first propose a formalism to describe the action of FA shimming pulses. Although a generalization would be possible for large flip angles, this theory for now is limited to the small tip angle regime where there exists a linear relationship between the RF excitation and the spin's tip angle [19].

Matrix
Representation of the Spatial Distribution the Transmit RF Field. As the control field distribution is typically obtained by MR acquisition, it is common in practice that ω 1 is not defined at each voxel of the underlying MR acquisition (existence of regions where the reconstruction of B + 1 is not possible). We thus define by R the set of voxel positions r where a given realization of ω 1 (r) is known (R is subject and session dependent). In the following, we shall use the following matrix representation of the spatial distribution of the control field: where M � |R| is the number of positions contained in R.
When there is no ambiguity about the choice of the region of interest, we shall simply use the notation Ω 1 .

Flip Angle Shimming
Pulse. Let u(t) ∈ C 1×N c and k(t) ∈ R 1×3 be the RF and k-space trajectory waveforms of a given FA shimming pulse. Let T be the duration of the RF and magnetic field gradient (MFG) waveforms and ω 0 (rad/ s) be the resonance offset at location r due to the static field nonuniformity. Let us use below the notation P(r, t) � u(t)e ik(t)·r . In the small tip angle approximation, the effective (time-averaged) control field created by this RF and magnetic field gradient (MFG) pulse pair is given by [8] where Note that the flip angle created by the pulse P is simply ϕ(P, r) � T|ω 1 (P, r)|. As a result, the mean, μ ω 1 , and variance, σ 2 ω 1 , of the effective control field ω 1 are given by μ ω 1 (r) � P(r)μ(r), In this model, the effect of off-resonances is taken into account but is not modelled as a random process ω 0 (r). Equation (4) hence can be implemented as such if the static field offset is mapped in a subject-specific manner. Alternatively, one can assume ω 0 � 0, i.e., that the static field offset has a negligible influence on pulse performance, which is acceptable for short (Tω 0 ≪ 1) pulses.
In the following, we present a method to compute universal FA shimming pulses which is entirely based on a second-order statistical model for the control field. Along with this statistically driven framework, we provide a method to evaluate the performance of a universal FA shimming pulse.

Statistical Robust Design of Universal Pulses.
Let P be a candidate FA shimming pulse to create a unitary effective control field ω 1 . Let us denote by the discrepancy between the unit target and |ω 1 (P, r)|, and by Π P (Δ, r) the probability that δ(P, r) exceeds a certain error level Δ > 0. To optimize P, we can exploit here a statistical robust design (SRD) approach [20]. Let us introduce a new design parameter 0 ≤ θ ≤ 1, called SRD tolerance, and define the error level Δ(P, θ) as where R design denotes the region of interest for the design of the pulse (any position outside r being therefore ignored in this metric). e SRD problem may then be posed as the minimization problem: Denoting now by F δ(P,r) the cumulative distribution function of δ(P, r), the condition Π P (Δ, r) ≤ θ gives Hence, the problem described by equation (8) takes the equivalent form: where F − 1 δ(P,r) denotes the inverse function of F δ(P,r) . A useful simplification of this optimization problem is then obtained by replacing the maximum (L ∞ norm) by the L q norm for a certain choice of q ≥ 1. is leads to the following new pulse optimization problem: If, now, for a given pulse candidate P, the following two conditions are satisfied for all positions r in R design : var ω 1 (P, r) ≪ ω 1 (P, r) 2 , then we can show (see Appendix A) that the pulse optimization problem (11) can be further simplified to take the form: where the dependence on the second-order B + 1 statistics can be made explicit: Retrospectively, the link with the database-driven objective function (simultaneous design across all subjects of the database of B + 1 maps), denoted here by f UP , can be explained as follows. e latter can be defined as the mean across the design database of the mean squares error of the effective control field: Now, if we assimilate sample mean and statistical mean, we obtain which leads us to the same expression as the one obtained in equation (13) for q � 2. is shows the near equivalence (note however that the objective function defined in reference [15] is proportional to the average (L 1 norm) of the FA-NRMSE rather than the L 2 norm and that the brain region for the evaluation of the FA-NRMSE was subjectspecific rather than fixed as in equation (15)) of the databasedriven and the statistical approaches when sample mean and statistical mean can be assimilated.
is is acceptable in particular when the design database is large enough to be representative of the whole targeted population.

Derivation of a Lower Bound for the Coefficient of Variation of the Flip Angle.
Let α(r) be a solution to the problem of minimizing the variance of the control field at position r among the RF shim solutions that return on average, up to a phase term, a unitary control field: In Appendix B, an explicit description of α(r) is obtained under the condition that e expression for α(r) then reads and the variance of the value of the objective at the optimum is given by Furthermore, it is shown in Appendix C that the locally optimal solution α(r) sets a performance limit for the coefficient of variation (CV) of the effective control field: which reads Hence, under the condition β ≫ 1, the lower bound for the CV of the effective control field, which we may also call ultimate FA-CV in the linear STA regime, and below denoted by g ∞ , reads e condition β ≫ 1, on the other hand, has a clear physical interpretation: the mean control field magnitude μ dominates its standard deviation. e validity of this condition in practice will be tested subsequently using experimental data.
It is also of interest to define the ultimate FA-CV-in the sense of equation (23)-for the circularly polarized (CP) mode of the TX array. In this mode, a fixed relationship between the RF signals applied to the TX array ports is imposed in order to reproduce the behavior of a single port volume coil. e CP mode can be defined with the vector a ⌣ ∈ C 1×Nc which describes the way the channels of the TX array interfere. e mean and variance of the CP mode control field, (μ Application of equation (23) on the CP mode control field statistics gives where In what follows, we present an application of the proposed statistical model of the intersubject variability of the control field and their implication in terms of pulse design and flip angle shimming performance assessment.

Control Field Second-Order Statistics.
e data consisted of a collection of N � 40 B + 1 maps measured on different adult subjects (age � 40 ± 15 years) which were recruited as part of previous studies, described in [14,15,21]. e MRI data were obtained on a parallel transmit-enabled whole-body Siemens 7 T system (Siemens Healthineers, Erlangen, Germany) equipped with an SC72 whole-body gradient insert (200 mT/ m/ms and 100 mT/m maximum gradient slew rate and maximum gradient amplitude, respectively) and the 8Tx-32Rx Nova head coil (Nova Medical, Wilmington, MA, USA) for signal transmission and reception. is RF coil array possesses a CP mode for transmission in which the amplitude and phases at excitation ports 2 to 8 are imposed by the amplitude and phase applied to the first excitation port. is operation mode thus has only one degree of freedom (the complex amplitude applied to the first excitation port) and is designed to create the TX profile corresponding to a birdcage resonator. e parallel transmission system consisted of eight transmitters (1 kW peak power per channel) driven with the Siemens Tim Tx (Step 2. 3) system. e protocol used to reconstruct the quantitative B + 1 maps was acquired from a multislice interferometric turbo-FLASH acquisition [22,23]  e control field second-order statistics (μ(r), C(r)) of the TX array head coil were computed from the TX field maps of the first (chronology of data acquisition) 35 subjects of the database using the mean and covariance estimators, as described in the eory section. For the estimation of the covariance matrix, the Ledoit-Wolf correction was used due to the limited number of subjects as compared to the dimension of the covariance matrix (64).

Coefficient of Variation of the Flip Angle.
e lower bounds for the CV of the effective control field ω 1 in the linear STA regime for the CP and pTX modes were computed from the estimated second-order statistics of the control field (μ, C) using equation (23). is enabled us in particular to compare the ultimate FA shimming performance when exploiting full pTX capability (independent RF signals applied to each port of the TX array) with the FA shimming performance of the CP mode of the TX array.

Statistical Robust Design of k T -Point Pulses.
For that study, the spins were assumed to be at resonance ω 0 � 0 everywhere, i.e., we assumed that the off-resonance effect on excitation was negligible. Computed FA shimming pulses were 7-k T -point [11] pulses (7 rectangular subpulses), wherein the duration of the RF subpulses and MFG blips were set to 140 μs and 50 μs, respectively. In all optimizations, the target flip angle, ϕ target , was 10 deg. Furthermore, the RF peak power applied at each excitation port was limited to 500 W, the total RF energy (respectively, the RF energy per TX channel) of the pulse was assigned to not exceed 60 mJ (respectively, 12 mJ), and the k-space displacement in x, y, and z created by each MFG blip was imposed to not exceed 50 rad/m.
A first FA shimming pulse, referred to as SRD-UP, was designed from the second-order control field statistics of the Nova TX array head coil (estimated in Section 3.1) by solving numerically the problem (same as equation (13) but nonunitary target, T duration of the pulse): In the above expression, R design was defined as the set of positions that belonged to the brain region in at least half of the subjects of the database, i.e., where R (n) brain denotes the mask of the brain for the n-th subject of the database and 1 R (n) brain (·) is the indicator function of this region.
For comparison with SRD-UP, the solution (same parameterization and constraints as for the design of SRD-UP) of the database-driven UP objective was also computed. e database-driven UP objective was defined as follows: is pulse is referred to as UP. Finally, a standard hard pulse where the transmit array is used in CP mode was defined. e amplitude of this pulse thus was set to return on average on the design database the target flip angle.
We then studied the spatial distribution of the mean and CV of the flip angle for UP, SRD-UP, and CP and compared it with the ultimate FA shimming performance g ∞ . For convenience, this performance is attributed to a virtual pulse P which would satisfy (see equations (4) and (19)) TP(r) � ϕ target α(r). In the following, this pulse is referred to as UP ∞ . Note that α(r) being a local solution, no RF and MFG waveforms can satisfy strictly this relationship unless the k-space trajectory defined by this pulse offers a complete coverage of the transmit k-space trajectory. is would lead in practice to a prohibitively long FA shimming pulse.
Using the five last subjects of the cohort (the 35 first constituting the design database), the flip angle normalized root mean squares error (FA-NRMSE) of UP, SRD-UP, and CP (P one of these pulses, ϕ Bloch (P, r), the FA created by pulse P at location r, 36 ≤ n ≤ 40): was computed on each subject. For the FA-NRMSE computation, the FA maps, ϕ Bloch (P, ·), were obtained by numerical integration of Bloch's equations (whereby T1 and T2 relaxations were neglected) rather than using the linear small tip angle approximation (equation (3) For every subject of this test database, a subject-tailored k T -point pulse (minimization of the n th term of the sum in equation (28)), possessing the same parameterization as the above universal pulses (SRD-UP and UP), was also designed to minimize the FA-NRMSE. ese solutions are referred to as ST 36 , . . ., ST 40 .
All pulse optimizations were performed using the second-order active-set algorithm, available in the Optimization Toolbox of MATLAB R2019a ( e Mathworks, Natick, MA), which allows solving constrained optimization problems numerically. e active-set algorithm used the analytical expressions for the Jacobian of the objective function and the RF power, RF energy, and MFG constraints (Hoyos-Idrobo et al. [24]). e number of iterations was limited to 2 × 10 3 . e algorithm however could stop before exceeding this number of iterations provided it reached a local minimum within default tolerance. Figure 1, the secondorder statistics of the control field for the Nova TX array are displayed for one axial slice through the brain. e mean control field μ thus gives rise to eight axial images, each one Concepts in Magnetic Resonance Part A corresponding to one TX channel. We also provide the mean control field of the CP mode, μ ⌣ � a ⌣ μ. e covariance matrix is composed of the control field variance of every TX channel (eight diagonal images) as well as the covariance between channels (off-diagonal images). Figure 2 displays sagittal, coronal, and axial images of the lower bound for the CV of the effective control field (ultimate FA shimming performance) in CP mode and the pTX operation mode.

Ultimate FA Shimming Performance.
Here, it can be seen that g ∞ < g ⌣ ∞ , i.e., that the pTX drive mode allows for a significant decrease in the CV of the effective control field as compared to the CP mode. In other words, for every voxel in the brain, there exists a pTX excitation mode that yields less intersubject variability than the CP mode. Quantitative comparisons of these performance boundaries are also provided in the form of cumulative histograms of g ∞ and g ⌣ ∞ in Figure 3. Finally, histograms of the distributions of the ratio of CVs g ⌣ ∞ /g ∞ are provided in Figure 4. From this plot, it can be seen in particular that the CV of the effective control field ω 1 can be reduced by a factor greater than two in some locations in the brain with pTX as compared to the CP mode.

FA Shimming Performance Simulation.
e FA shimming performance analysis of UP ∞ , UP, SRD-UP, and for the standard hard pulse in the CP mode is shown in Figures 5  and 6. For UP and SRD-UP, in the brain (in the sense of (30)), the mean flip angle ϕ(r) (statistical mean) is very close to the target (10deg) as it remains in the interval [9,11]deg. By definition of α (see equation (19)), for the virtual pulse UP ∞ , the mean flip angle is exactly equal to the target. As it can be expected, for the CP mode, the mean flip angle can deviate very strongly from the 10deg target as it reaches for instance in the center >20deg. On the other hand, the median CV of the flip angle amounts to 12% for UP and CP and 10% for SRD-UP and is ultimately <8% for UP ∞ . In the "worst location" which was found to be the cerebellum (see Figure 5), the CV of the flip angle does not exceed 15% with UP ∞ but can exceed 20% with "real" pulses. e result of the FA-NRMSE analysis made separately on the design and test subjects is presented in Figure 7. e results obtained on the test subjects are compliant with the FA-NRMSE statistics obtained on the design database in the sense that no significant offset (unpaired Wilcoxon ranksum test) was found between the FA-NRMSE obtained in the design and test databases. On the test subjects, the performance was very comparable between UP ∞ , UP, and SRD-UP, with a median FA-NRMSE of 5 to 7%. e subjecttailored pulses returned FA-NRMSEs <2%.

Storage and Computation Time Efficiency. Denoting by
N v the average number of voxels in the brain, the storage volume for the first-and second-order statistics of the control field and for a database of N subjects are (N c × (N c + 3) × N v )/2 and N × N v × N c , respectively. us, for N ≥ (N c + 3)/2 (6 in our case), the statistical representation becomes more parsimonious. e computation time for one evaluation of the objective function (and the Jacobian) was 2.13 s and 0.26 s for UP and SRD-UP, respectively. As these workloads scaled linearly with the number of voxels (4.22 × 10 5 in total for the 35 subjects of the database and 1.1 × 10 4 for R design (see equation (27))), we found that the workloads per voxel were 5 μs for UP and 17.4 μs for the computation of UP and SRD-UP, respectively. As a result, but keeping in mind that these numbers may depend on the pulse parameterization (e.g., number of k T -points), the design of SRD-UP is computationally more efficient than the UP design when the size of the design database is greater than 17.4/5 ≃ 3.5.

Discussion
In this work, we have shown that the statistical robust and the UP design methods are theoretically nearly equivalent when the size of the database is sufficiently large to embrace the variability of the resonant TX RF field across the population of interest. Despite its apparent complexity, the statistical robust design algorithm defined by equation (13) has a rather straightforward implementation as the general term of the summation across voxels depends explicitly on the mean vector and covariance matrix of the TX RF field (see equation (14)). As compared to the database-driven UP design methodology (Gras et al. [14]), an important difference is the involvement of the covariance matrix C(r) which has not equivalent in the original formulation. To our knowledge, the design of UPs using the SRD formalism, originally applied by Alotto and colleagues for the design of electromagnetic structures, has not been proposed so far for RF pulse design and could represent some advantages. First, as the SRD computation time does not depend on the size (number or subjects) of the database, this approach enables integration of potentially very large databases, which is a priori necessary to ensure optimal UP performance. Secondly, equipped with this model, one can likewise compute probabilities for obtaining a certain deviation from a target flip angle. Doing such estimations using a database of TX field maps would require unrealistically large number of subjects.
e theory presented in this work defines the performance of a FA shimming pulse as the coefficient of variation (std/mean) of the flip angle. Furthermore, in the linear STA regime, it is shown that this CV cannot be lower than a certain value g ∞ which depends on the local second-order statistics of the control field. In the process of validating-and ultimately certifying-the use of pTX in a clinical setup, this performance measure may prove useful. One concrete application of this concept is the comparison of the CP mode ultimate performance (g ⌣ ∞ ) with the pTX performance (g ∞ ): for the Nova TX array, we found that the availability of eight independent TX channels to locally tune the drive mode so as to minimize the variance of the effective control allowed reducing the CV of the flip angle by a factor of ≃ 2, while achieving close to the target on average (see Figure 4). Interestingly, this theoretical analysis sheds some light onto the concept of UPs. Indeed, with statistical arguments, one can take advantage of the TX array architecture to mitigate both the spatial heterogeneity and the intersubject variability of the transmit RF field in the brain. Finally, the coefficient of variation metric for UP ∞ (see equation (23)), being characteristic of the TX array architecture, can potentially serve as a guide for coil-design purposes. e present study was performed with the B + 1 data collected on 40 subjects, from which a subgroup of 35 subjects was used to compute statistics. Given the high dimension of the covariance matrix of the control field for an eightfold TX array system, it was necessary to use the Ledoit-Wolf correction for the estimation of the covariance matrix of the control field (to obtain positive definite estimates). However, we can still expect a bias in the estimation of the covariance matrix. As a result, it is possible that the CV analysis reported here (lower bound and CV maps for the different k T -point pulses computed in this study) can  Figure 2: Sagittal, coronal, and axial views of the lower bound for CV of the effective control field (ω 1 ) for the CP drive mode (left), g ⌣ ∞ , and the full pTX drive mode (right), g ∞ . As compared to the CP drive mode, FA shimming performances, measured as the CV of the effective control field on each location, are clearly superior using pTX (i.e., 8 degrees of freedom as opposed to 1 for the CP mode). It also shows that in CP mode, the greatest intersubject variability occurs in the cerebellum and in the temporal lobes. deviate from the analysis that one would obtain with a larger database to estimate the control field second-order statistics. e theory developed in this study is based on the assumption that the statistical distribution of the control field vector ω 1 (r) is Gaussian. is assumption is very convenient from the mathematical point of view as this property is inherited by the effective control field ω 1 (r) in the linear STA regime. With a few tens of subjects, this hypothesis has not yet been verified quantitatively, but given the complexity and the multiplicity of effects leading to the observed variability, the Gaussian model is a reasonable assumption a priori. Another hypothesis made in this study is the fact that β � μ * C − 1 μ is much greater than one. Here we notice that this is equivalent to the assumption that g ∞ is much smaller than one. In Figure 2, g ∞ was seen not to exceed 0.2, confirming that this hypothesis was approximately valid over the whole brain.
Extending this work to large flip angle pulse is one interesting prospect. However, in this case, the linear relationship between ω 1 and ω 1 exploited to derive a simple g ∞ (CP)/g ∞ (pTX) Figure 4: Cumulative histogram of the ratio g ⌣ ∞ /g ∞ , whereby g ∞ and g ⌣ ∞ are the theoretical lower bounds for the effective control's CV in the pTX (8 degrees of freedom) and the CP (1 degree of freedom) modes. e more "agile" pTX drive mode can adapt to the known control field statistics and thereby allows reducing significantly the lower bound for the CV of the effective control field. A reduction factor of at least 1.5 is obtained for 80% of the voxels of the brain. In some areas, the CV reduction factor reaches a factor of 2 or more.     expression of the SRD objective does not hold. Here a pragmatic approach would be to replace exact calculation of the first and second moments of the flip angle (seen as a random variable) by numerical approximations thereof. A possible direction to perform such calculation is the use of the so-called unscented transform method [25] which seems here particularly adapted [20].
An important difference with the pulse design method reported in [13] is the deterministic nature of the off-resonance term ω 0 in the current statistical design. In this study, we set ω 0 � 0 for the design of UPs using the SRD formalism, i.e., we neglected off-resonance effects in the performance assessment. While acceptable for short excitation pulses, this can be a significant source of error in practice when the pulse duration becomes large [26]. Extending the statistical robust design approach to account for a random off-resonance term seems however feasible with some additional hypotheses. Hence, if it is assumed that ω 0 (r) is a centered Gaussian random variable independent of ω 1 , it remains possible in equation (3) to compute the mean deviation of the effective control field with respect to the effective control field at resonance. is quantity takes the form: where ω (0) 1 denotes the effective control field at resonance, and where the expectation operator · 〈 〉 covers ω 1 as well as ω 0 statistics. A possible direction is to exploit this quantity as an additional objective to improve the robustness against off-resonances while taking into account the regional dependence of the variance of the resonance frequency offset.
In this work, for simplicity, the statistical robust pulse design implementation did not include specific absorption rate (SAR) constraints. However, as it was done in the past for the design of UPs, there would be no obstacle in adding explicit local SAR constraint to the statistical robust design procedure in the form of virtual observation points [27]. Interestingly also, the RF electric field distribution (E-field) at the origin of the SAR is another imperfectly controlled physical quantity affecting RF safety management procedure, in particular at UHF and with transmit array architectures [28]. Like the B + 1 field, the E-field could advantageously be treated as well as a stochastic process, thereby allowing for a probabilistic model of the local 10 g SAR. Naturally, the dual problem consisting of minimizing the SAR under performance (flip angle error) constraints [29] could benefit as well from such a unified statistical model. One obstacle for the estimation of the E-field statistics yet remains the absence of direct method to measure the E-field in vivo and thus the necessity to resort to electromagnetic simulations on a large amount of head models to give access to them. Furthermore, due to the interweaved nature of the RF electric and magnetic fields, the values of the Eand B + 1 -field, seen here as random processes, are in fact necessarily correlated.

Conclusions
We have presented in this work a statistical robust approach to design calibration-free pTX universal pulses for brain imaging at ultrahigh field. Experimentally, this method has been shown to provide comparable FA-NRMSE results as the conventional database-driven UP design, but presents a significant advantage in terms of computational efficiency when the size of the design database exceeds 4. Furthermore, in contrast with the database-driven design of UP, the computation time of the statistical robust UP does not depend on the size of this database and consequently enables integration of a very large design database. is is an important criterion to ensure optimal UP performance. Furthermore, the theoretical lower bound for the coefficient of variation of the flip angle derived in this work, g ∞ , appears as a valuable analytical performance limit to assess the FA shimming performance of a given pulse. Finally, this metric can be useful to characterize, and possibly also optimize, the performance of TX array architectures.

A. Statistical Robust Pulse Design
In this section, we focus on the relaxed SRD problem (infinite norm replace by the L q norm) given by equation (11). Let us study the first-and second-order statistics of δ(P, r).
From ω 1 ∼ N(μ, C), we have: ω 1 � Pω 1 ∼ N(Pμ, PCP * ) and |ω 1 | Rician. Now, assuming that for all voxel positions, the mean effective control largely exceeds its standard deviation, i.e., we can assume that ω 1 owns a Gaussian distribution with the following mean and variance: It follows from condition (A.1) that the distribution of δ(P, r) � ‖|ω 1 | − 1 | is the so-called folded normal distribution and that the mean, μ δ , and the variance, σ 2 δ , of δ(P, r) are therefore given by var ω 1 e − ξ 2 /2 , where ξ is defined as follows: and where Φ denotes the cumulative density function of the standard normal distribution N(0, 1). In Figure 8, we have plotted for three values of the SRD design variable θ the quantity: as a function of ξ. e graph of κ θ as a function of ξ is very flat in ξ � 0 and, for |ξ| < 1 and θ ≥ 0.01, we have (see Figure 8) As a result, for a candidate pulse P for which condition (A.6) is satisfied for all r ∈ R design , then the objective in equation (11) simplifies to Noting that (i) the leading term in equation (A.7) is independent of r and (ii) that the pulse design problem thus simplifies to min P r∈R design We note that the condition given by equation (A.6) supposes that the candidate FA shimming pulse P returns an effective control field that is "almost" uniformly unbiased ( |ω 1 | ≃ 1 everywhere). erefore, it would not be realistic to assume that this condition is fulfilled for all P. But, for the above simplification to be made, one only need to assume that this condition is fulfilled in the neighborhood of the solution of (A.9).

B. Lower Bound for the Coefficient of Variation of the Effective Control
Below we first solve the phase-constrained problem and then the phase-unconstrained problem. While the first problem can be solved analytically in general, the resolution of the second problem-the one that we focus on in this study-is conducted under the assumption that β ≫ 1. Using the equality constraint, we have e Lagrangian of the above problem reads L � aCa * + λ(aμ − 1).

(B.3)
If a is solution, then where ρ(a) � |aω 1 |. Assuming Gaussian statistics for ω 1 (ω 1 ∼ N(μ, C)), then aω ∼ N(aμ, aCa * ) and so ρ(a) is Rician distributed with mean (L 1/2 (x) � e x/2 ((1 − x)I 0 (− (x/2)) − xI 0 (− (x/2))) where I 0 and I 1 are the modified Bessel functions of the first kind with parameter α � 0 and α � 1, respectively): Given now that aμμ * is colinear to μ * , a necessary condition to have ∇ a L � 0 is that aC is colinear with μ * , and so, that there exists a scalar t such that a � tμ * C − 1 . From equation (B.16), we obtain (see equation (18) As a result, the solution to the phase-unconstrained problem is given (up to a uniform phase term across TX channels) by Taking the above expression for a in equation (B.11) (the condition under which the statistical distribution for ρ(a) can be assumed Gaussian), and observing that |αμ| 2 � (β 2 /t 2 ) and αC − 1 α * � (β 2 /t), we obtain that the condition to be satisfied so that the distribution for ρ(α) is Gaussian is β ≫ 1. Under this condition, we observe that the optimum of the phase-unconstrained problem converges to the optimum of the phase-constrained problem.

C. Locally Optimal FA Shimming Configuration
Given a FA shimming pulse P � u ⊙ k, and a given position r 0 , let us define Q r 0 (r, t) as follows: Q r 0 (r, t) � P(r, t) P r 0 , t ω 1 r 0 .

(C.5)
To complete the proof, we note that the function z ↦ |z| is a convex function and so that for any random complex random variable Z, |Z| 〈 〉 ≥ | Z 〈 〉|. Applying this to Z � P(r 0 , t)ω 1 (r 0 ) yields var P r 0 , t ω 1 r 0 P r 0 , t ω 1 r 0 2 ≥ var α r 0 ω 1 r 0 . (C. 6) Data Availability e datasets generated and/or analyzed during the current study and the pulse design source codes are available from the corresponding author on reasonable request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.