# Rotating vector model for magnetars

H. Tong1, P. F. Wang2, H. G. Wang1, Z. Yan3
1School of Physics and Materials Science, Guangzhou University, Guangzhou
510006, China
2National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun
Road, Chaoyang District, Beijing 100012, China
3Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, China E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The modification of the rotating vector model in the case of magnetars are
calculated. Magnetars may have twisted magnetic field compared with normal
pulsars. The polarization position angle of magnetars will change in the case
of a twisted magnetic field. For a twisted dipole field, we found that the
position angle will change both vertically and horizontally. During the
untwisting process of the magnetar magnetosphere, the modifications of the
position angle will evolve with time monotonously. This may explain the
evolution of the position angle in magnetar PSR J1622-4950 and XTE J1810-197.
The relation between the emission point and the line of sight will also
change. We suggest every magnetospheric models of magnetars also calculate the
corresponding changes of position angle in their models. Order of magnitude
estimation formula for doing this is given. This opens the possibility to
extract the magnetic field geometry of magnetars from their radio polarization
observations.

###### keywords:

stars: magnetars – pulsars: individual (PSR J1622-4950, XTE J1810-197) –
pulsars: general

††pubyear: 2015††pagerange: Rotating vector model for magnetars–C

##  1 Introduction

Magnetars are thought to be young and strongly magnetized neutron stars (Kaspi
& Beloborodov 2017). There may be a continuous transition from normal pulsars
to magnetars. Observational evidences include: (1) Five magnetars are observed
to have pulsed radio emissions (Camilo et al. 2006, 2007a; Levin et al. 2010;
Eatougth et al. 2013; Lower et al. 2020a). (2) Magnetar-like activities are
observed in two high magnetic field pulsars and one is also radio-loud
(Gavriil et al. 2008; Archibald et al. 2016). (3) Several low dipole magnetic
field magnetars are also discovered (Rea et al. 2010; Zhou et al. 2014) etc.
The burst and persistent X-ray emissions of magnetars are thought to be
powered by the neuton star’s strong magnetic field, i.e. the name “magnetar”
(Duncan & Thompson 1992). The reason why magnetar can have a huge magnetic
energy release rate may be that their magnetic field is twisted compared with
that of normal pulsars (Thompson et al. 2002; Beloborodov 2009; Paven et al.
2009; Tong 2019). The untwisting of the twisted magnetic field is responsible
for the outburst activities seen in magnetars (Coti Zeltati et al. 2018).

Up to now, all the evidences for a twisted magnetic field in magnetars are
indirect observations (Tiengo et al. 2013; Weng et al. 2015). Direct evidence
of a twisted magnetic field may be revealed by future X-ray polarimetry
observations. At the same time, for radio emitting magnetars, we already have
polarization observations which may tell us the magnetic field geometry of
magnetars (Camilo et al. 2007b; Kramer et al. 2007; Camilo et al. 2008; Levin
et al. 2012; Eatought et al. 2013; Lower et al. 2020b). However, when
analysing polarization position angle of the radio emitting magnetars, the
rotating vector model for a dipole field is always employed (e.g., Levin et
al. 2012; Lower et al. 2020b). The rotating vector model in the case of a
twisted magnetic field is not available at present. This hinders us from
extracting the magnetic field geometry information from magnetar radio
observations.

We will build a rotating vector model in the case of twisted dipole field for
magnetars. It is compared with the current magnetar radio observations. It is
shown that the polarization of magnetar radio emission may be modeled to some
degree by the twisted dipole field.

The structure of this paper is: the rotating vector model for a twisted dipole
field is calculated using simple geometry in Section 2. A more detailed
calculation is presented in Section 3. A twisted dipole field will also change
the relation between the emission point and the line of sight. This is shown
in Section 4. Comparison with the observations is presented in Section 5.
Discussion and conclusion are given in Section 6. In normal pulsars, it is a
common method to extract magnetic field geometry using the rotating vector
model. The corresponding calculations in the dipole case are given in the
Appendixes.

##  2 Calculations using simple geometry

The geometry of a rotating dipole is shown in figure 1. The angle between the
magnetic axis and the rotational axis is denoted as $\alpha$ (inclination
angle), the angle between the line of sight and the rotational axis is denoted
as $\zeta$ (viewing angle). The plane where the rotational axis and the
magnetic axis lies is denoted as $\phi=0$ (the meridian plane). The closest
approach between the line of sight and magnetic axis is denoted as $\beta$
(impact angle). From figure 1, it can seen that $\beta=\zeta-\alpha$. The
position angle $\angle MPR=\psi$ can be found as a function of the rotational
phase $\angle MRP=\phi$

$\tan\psi=\frac{\sin\alpha\sin\phi}{\cos\alpha\sin\zeta-\sin\alpha\cos\zeta\cos\phi}.$
(1)

The deduction process is given in Appendix A. Polarization position angle
curve for typical parameters is shown in figure 2. A typical inclination angle
($\alpha=45^{\circ}$) and a large impact angle ($\beta=20^{\circ}$) is chosen
in figure 2. Please note that the position angle defined in equation (1) and
the observationally reported position angle is different by a minus sign:
$\psi_{\rm obs}=-\psi$ (Everett & Weisberg 2001). This will account for the
position angle of PSR J1622-4950: it has a positive slope in figure 4 in Levin
et al. (2012), while the reported impact angle $\beta$ is negative in Levin et
al. (2012).

Figure 1: Geometry of a rotating dipole. The rotational axis and the magnetic
axis are shown, along with the line of sight. At the observational point $P$,
the magnetic field line makes a position angle $\psi$ with the rotational
axis. The red solid line is for a twisted dipole field. Figure 2: Polarization
position angle as a function of pulse phase. The black solid line is for a
dipole magnetic field, with $\alpha=45^{\circ}$, $\zeta=65^{\circ}$
($\beta=20^{\circ}$). The red dashed line is for a twisted dipole field, with
twist parameter $n=0.5$, the inclination angle and viewing angle is the same
as the dipole case.

For a twisted dipole field, the magnetic field will have a toroidal component
compared with the dipole case (Wolfson 1995; Thompson et al. 2002; Pavan et
al. 2009; Tong 2019). The toroidal component will result in modifications to
the position angle. For small twist, the magnetic field can be considered as a
dipole field (which is purely poloidal) plus the toroidal component
(Beloborodov 2009; Tong 2019). Therefore, the corresponding position angle may
also be considered as the dipole case plus some modification

$\psi=\psi(\rm dipole)+\Delta\psi_{\rm twist},$ (2)

where $\psi(\rm dipole)$ is the corresponding value given in equation (1), and
$\Delta\psi_{\rm twist}$ is due to the twist of magnetic field lines. For a
self-similar twisted dipole field with radial dependence $B(r)\propto
r^{-(2+n)}$, the twist angle will be the same for the same initial and final
magnetic colatitude. Here the parameter $n$ characterize the twist of the
magnetic field: $n=1$ corresponds to the dipole case, $n=0$ the split monopole
case, and $n\in(0,1)$ is the general twisted dipole case. When referring to
the magnetic field, we will mainly work in the magnetic frame, where the
magnetic axis is the spherical polar axis. For an initial magnetic colatitude
$x_{1}=\cos\theta_{1}$, final magnetic colatitude $x_{2}=\cos\theta_{2}$, the
twist of magnetic field is (Wolfson 1995)

$\Delta\phi=\frac{\lambda}{n}\int_{x_{2}}^{x_{1}}\frac{f^{1/n}(x)}{1-x^{2}}dx,$
(3)

where

$\lambda=\sqrt{(35/16)(1-n)},$ (4)

and $f(x)$ is the dimensionless flux function. About the twist direction of
the magnetic field, a positive $\lambda$ corresponds to eastward motion of
magnetic field lines in the southern hemisphere (Tong 2019). See figure 3 for
an example. The twist field lines are also shown in figure 2 in Wolfson (1995)
and figure 3 in Thompson et al. (2002). The modification to the position angle
is positive when the field line twist westward (as shown in figure 1, negative
$\lambda$). Otherwise, it will be negative (positive $\lambda$).

The maximum twist is defined as the twist along the magnetic field line from
the north pole to the sourth pole

$\displaystyle\Delta\phi_{\max}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{n}\int_{-1}^{1}\frac{f^{1/n}(x)}{1-x^{2}}dx$ (5)
$\displaystyle=$
$\displaystyle\frac{2\lambda}{n}\int_{0}^{1}\frac{f^{1/n}(x)}{1-x^{2}}dx$ (6)
$\displaystyle\approx$ $\displaystyle 2\lambda,$ (7)

where during the calculations, the small twist approximation is used $n\approx
1$, and $f(x)\approx 1-x^{2}$. This approximation is accurate for a wide
parameter space (Tong 2019). Therefore, there is a one to one corresponds
between the three paprameters: $n$, $\lambda$, and $\Delta\phi_{\rm max}$. For
the case shown in figure 1, the modification to the position angle is

$\displaystyle\Delta\psi_{\rm twist}$ $\displaystyle=$
$\displaystyle\frac{\lambda}{n}\int_{\theta_{\rm
obs}}^{1}\frac{f^{1/n}(x)}{1-x^{2}}dx$ (8) $\displaystyle\approx$
$\displaystyle\frac{\lambda}{n}(1-\cos\theta_{\rm obs})$ (9)
$\displaystyle\approx$ $\displaystyle\frac{\lambda}{2n}\sin^{2}\theta_{\rm
obs},$ (10)

where during the last step, the approximation of small magnetic colatitude is
used $1-\cos\theta_{\rm obs}\approx 1/2\sin^{2}\theta_{\rm obs}$. Adding the
correct sign of the twist angle, the final result is

$\displaystyle\Delta\psi_{\rm twist}$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{n}(1-\cos\theta_{\rm obs})$ (11)
$\displaystyle\approx$ $\displaystyle-\frac{\lambda}{2n}\sin^{2}\theta_{\rm
obs}.$ (12)

As discussed above, a positive $\lambda$ corresponds to a magnetic field twist
eastward. The position angle will be smaller. This corresponds to the negative
sign in the equation (12). In below, equation (12) is confirmed using more
detailed calculations.

From equation (12), it can be seen that the modification to the position angle

$\Delta\psi_{\rm twist}\propto\lambda\propto\sqrt{1-n}\propto e^{-t/2\tau},$
(13)

where during the last step the evolution of the parameter $n$ is adopted from
the toy model for magnetar outburst (Tong & Huang 2020). The parameter $\tau$
is the exponential decay timescale of the X-ray flux or hot spot. Therefore,
the modification to the position angle, especially the constants of the
position angle $\psi_{0}$, will evolve with time monotonously. This is
consistent with the polarization observation of the radio emitting magnetar
PSR J1622$-$4950 (Levin et al. 2012). Quantitative comparison will be given in
below.

Figure 3: Twisted magnetic field lines. The black solid line is for a dipole
field. The red solid line is for a twisted dipole field. It twists westward in
the southern hemisphere. The red dashed line twists eastward in the southern
hemisphere.

##  3 More detailed calculations

The above calculations may seem too simple minded. A more detailed calculation
can be done following the procedure of Hibschman & Arons (2001). The details
and corresponding calculations for the dipole case are given in Appendix B.
Considering the rotation of the neutron star (Blaskiewicz et al. 1991), or the
magnetospheric current (Hibschman & Arons 2001), the change of position angle
may be obtained using perturbations. These two effects are both proportional
to $r/R_{\rm L}$, where $r$ is the radio emission height, and $R_{\rm L}$ is
the light cylinder radius. For a twisted magnetic field, the result will be
qualitatively different. In the case of a self-similar twisted dipole field,
the modification to the position angle is independent of the emission height.
This is because all the field component (poloidal and toroidal) have the same
radial dependence $\propto r^{-(2+n)}$ (Wolfson 1995).

For an axisymmetric force-free magnetic field, the magnetic field obeys the
Grad-Shafranov equation. In the magnetic frame, using spherical coordinate,
the magnetic field is (Wolfson 1995)111A normalization factor for the magnetic
field omitted. It can be normalized to the polar or equatorial surface
magnetic field. Only the geometry of the magnetic field is included.

${\bf B}=\frac{1}{r\sin\theta}[\frac{1}{r}\frac{\partial
A}{\partial\theta}\hat{r}-\frac{\partial A}{\partial
r}\hat{\theta}+F(A)\hat{\phi}].$ (14)

For the self-similar case, the two functions are: $F(A)=\lambda A^{1+1/n}$,
$A(x)=r^{-n}f(x)$ (where $x=\cos\theta$). The toroidal magnetic field
component may be considered as perturbation to the dipole magnetic field
(Beloborodov 2009; Tong 2019). Therefore, the perturbation to the tangent
vector is

$\displaystyle{\bf t}_{1}$ $\displaystyle=$ $\displaystyle\frac{B_{\phi}}{|\bf
B|}\hat{\phi}\approx\frac{B_{\phi}}{|{\bf B}_{\rm
p}|}\hat{\phi}=\frac{B_{\phi}}{(B_{r}^{2}+B_{\theta}^{2})^{1/2}}\hat{\phi}$
(15) $\displaystyle=$ $\displaystyle\frac{F(A)}{[(1/r\,\partial
A/\partial\theta)^{2}+(\partial A/\partial r)^{2}]^{1/2}}\hat{\phi}$ (16)
$\displaystyle=$ $\displaystyle\lambda\frac{f^{1+1/n}(x)}{[(\partial
f/\partial\theta)^{2}+n^{2}f^{2}(x)]^{1/2}}\hat{\phi}.$ (17)

In order to an meaningful analytical result, further simplifications are
needed. At the same time, it is required that when the parameter $n$ is back
to one, the result should be back to the dipole case. Taken the flux function
as the dipole case $f(x)\approx 1-x^{2}=\sin^{2}\theta$, the parameter $n$ in
the power law index of $f(x)$ set to be $n=1$, while keeping the expression of
$\lambda$ unchanged, equation (17) can be further simplified

$\displaystyle{\bf t}_{1}$ $\displaystyle=$
$\displaystyle\lambda\frac{\sin^{3}\theta}{[4\cos^{2}\theta+\sin^{2}\theta]^{1/2}}\hat{\phi}$
(18) $\displaystyle\approx$
$\displaystyle\frac{\lambda}{2}\sin^{3}\theta\hat{\phi},$ (19)

where during the last step the approximation of small $\theta$ is used.

Following the procedure outline in Appendix B, it is straight forward (though
tedius) to calculate the changes in the position angle due to a twisted dipole
field

$\displaystyle\Delta\psi_{\rm twist}$ $\displaystyle=$
$\displaystyle-2\lambda\sin^{2}\theta$ (20) $\displaystyle\approx$
$\displaystyle-\frac{8}{9}\lambda\sin^{2}\theta_{\rm obs},$ (21)

where during the last step the relation between the line of sight and the
emission point are used: $\theta\approx 2/3\,\theta_{\rm obs}$, $\theta$ is
the colatitude of the emission point, $\theta_{\rm obs}$ is the colatitude of
the line of sight, both in the magnetic frame. $\theta_{\rm obs}$ can be found
as a function of the pulsar rotational phase, see equation (31). For positive
$\lambda$, the magnetic field will twist eastward. This will result in a
smaller position angle. This corresponds to the minus sign in equation (21).
The result in equation (21) is consistent with that in equation (12). The
observational meaning of equation (21) is that a twisted dipole field will
change the polarization position angle in the vertical direction. Figure 2
shows the modified position angle for a twisted dipole field, for twist
parameter $n=0.5$. A relatively large impact parameter is chosen. A large
impact angle will result in larger deviation from the dipole case. For a self-
similar twisted dipole field, there is one additional parameter $n$ (or
$\lambda$, or $\Delta\phi_{\rm max}$, there is a one to one correspondence
between these three parameters) compared with the dipole case. This is the
merit of assuming a self-similar twisted dipole field.

The relation between the emission point and the line of sight ($\theta\approx
2/3\,\theta_{\rm obs}$) is for the dipole case (Hibschman & Arons 2001). For a
twisted dipole field, the relation between the emission point and the line of
sight is only slightly different, as detailed in below.

##  4 Relation between emission point and line of sight

The relation between the emission point and the line of sight for a dipole
field is shown in Appendix C. For a twisted dipole field, the tangent vector
can be get from equation (14)

${\bf\hat{t}}=\left(\frac{1}{\sin\theta}\frac{\partial
f}{\partial\theta}\hat{r}+n\frac{1}{\sin\theta}f(x)\hat{\theta}+\lambda\frac{1}{\sin\theta}f^{1+1/n}(x)\hat{\phi}\right)\frac{1}{N},$
(22)

where $N$ is the normalization factor. For the small twist case, the
expression for the tangent vector can be further simplified

${\bf\hat{t}}=\left(2\cos\theta\hat{r}+\sin\theta\hat{\theta}+\lambda\sin^{3}\theta\hat{\phi}\right)\frac{1}{N}.$
(23)

In the small twist case, the normalization factor is
$N=(4\cos^{2}\theta+\sin^{2}\theta+\lambda^{2}\sin^{6}\theta)^{1/2}\approx(4-3\sin^{2}\theta)^{1/2}$.
In comparison with the normalization factor defined in Appendix C, here the
normalization factor is $N\approx 2N_{1}$. The main deviation from the dipole
case is a toroidal magnetic field component. The tangent vector can be
expressed in Cartesian coordinate. By equaling the tangent vector and the unit
vector along the line of sight, it is found that

$\displaystyle\frac{1}{N}(3\sin\theta\cos\theta\cos\phi-\lambda\sin^{3}\theta\sin\phi)=\sin\theta_{\rm
obs}\cos\phi_{\rm obs}$ (24)
$\displaystyle\frac{1}{N}(3\sin\theta\cos\theta\sin\phi+\lambda\sin^{3}\theta\cos\phi)=\sin\theta_{\rm
obs}\sin\phi_{\rm obs}$ (25)
$\displaystyle\frac{1}{N}(2\cos^{2}\theta-\sin^{2}\theta)=\cos\theta_{\rm
obs}.$ (26)

From equation (26), it can be inferred that the changes in $\theta$ is only
quantitative. The relation between $\theta$ and $\theta_{\rm obs}$ may be
taken as the same as the dipole case $\theta=(2/3)\theta_{\rm obs}$. The major
changes is for the azimuthal angle: $\phi=\phi_{\rm obs}+\delta\phi$. By
substituting into equation (25), keeping only the first order term, it is
found that

$\delta\phi=-\frac{4}{27}\lambda\sin^{2}\theta_{\rm obs}.$ (27)

For a positive $\lambda$, the field will twist eastward. The azimuthal angle
$\phi$ will be smaller than $\phi_{\rm obs}$. Therefore, the minus sign in the
above equation is correct. Equation (24) will give the same result. The
observational meaning of equation (27) is that a twisted dipole field will
change the polarization position angle in the horizontal direction.

By comparing equation (21) and equation (27) it can be seen that a twisted
dipole field will change the polarization position angle in both the vertical
and the horizontal direction. Furthermore, the changes in the vertical and the
horizontal direction have the same dependence on the twist parameter $\lambda$
and the colatitude of the line of sight in the magnetic frame $\theta_{\rm
obs}$. Quantitatively, the change in the horizontal direction is only $1/6$ of
that in the vertical direction $\delta\phi=1/6\ \Delta\psi_{\rm twist}$.
Therefore, the change of position angle in the vertical direction may be
easier to observe.

##  5 Comparison with the observations

The rotating vector model is commonly used in the case of normal pulsars to
obtain the information of the magnetic field geometry (Lyne & Manchester 1988;
Rankin 1993; Manchester et al. 1998; Johnston & Weisberg 2006; Johnston &
Kramer 2019). The position angle of some pulsars may be too complicated to be
modeled by the rotating vector model.

In the case of magnetars, previous polarization observations of five radio
emitting sources may already tell us some aspect of the magnetic field
geometry of magnetars. In our opinion, current observations already give us
some hints that the magnetic field of magnetars may be twisted. And the
rotating vector model for a twisted dipole field may explain the observations
to some degree.

  1. 1.

Changing $\psi_{0}$. For the third radio emitting magnetar PSR J1622-4950, its
position angle changes monotonously with time (figure 4 and table 1 in Levin
et al. 2012). For a twisted dipole field, the position angle will change in
the vertical direction (see equation (21)). This is also true for the position
angle in the meridian plane $\psi_{0}$. During the untwisting process, the
magnetic field will relax to that of the dipole case (Thompson et al. 2002;
Tong 2019). The parameter $n$ will approach 1. The parameter $\lambda$ will
decrease monotonously with time (equation (4)). According to the toy model for
magnetar outburst (Tong & Huang 2020), the parameter $\lambda$ will decrease
with time exponentially. The corresponding position angle, also for
$\psi_{0}$, will decrease with time exponentially (equation(13)). This is
consistent with the observations of PSR J1622-4950 (Levin et al. 2012). Figure
4 shows quantitatively that the $\psi_{0}$ of PSR J1622-4950 may decrease in
an exponential form222Mathematically, the observational data may also be
fitted using other forms, e.g. a linear decay form etc. According to our
physical modeling, we prefer an exponential decay form.. A decreasing
$\psi_{0}$ is also seen in the fifth radio emitting magnetar Swift
J1818.0-1607 (Lower et al. 2020b, table 3 there). Again, this indicates an
untwisting magnetic field in the putative magnetar.

From the typical exponential decay form in figure 4, the decay time scale is
about $90$ days. This is smaller than the X-ray flux decay time scale (360
days, Anderson et al. 2012). But they can be considered as the same order of
magnitude. The typical decay amplitude is about $25^{\circ}$. According to
equation (21), the changes of position angle at the meridian plane is about:
$|\Delta\psi_{0}|=(8/9)\lambda\sin^{2}\beta$, where $\beta$ is the impact
angle. For typical twist parameter $n=0.5$, the required impact angle for PSR
J1622-4950 in the twisted dipole scenario is about: $|\beta|=40^{\circ}$.
Observationally, the range of impact angle for PSR J1622-4950 is:
$\beta\in[-25^{\circ},-8^{\circ}]$ (Levin et al. 2012). The theoretical value
of the impact angle $\beta$ is slightly larger than the observational values.
But again, they can be considered as roughly consistent. In the presence of
multipole field, the modification to the twist angle may be larger.

  2. 2.

Changing slope. For the radio emitting magnetar PSR J1622-4950, the slope of
the position angle is also changing with time (Section 5.3 in Levin et al.
2012). A changing slope of the position angle is also seen in the first radio
emitting magnetar XTE J1810-197 (Section 3.4 in Kramer et al. 2007). During
the untwisting process of the magnetar’s magnetic field (Thompson et al. 2002;
Tong 2019), the parameter $n$, $\lambda$ etc are all evolve with time
monotonously. Then the change in position angle due to the magnetic field
twist will also evolve with time. This may explain the changing slope of
position angle in XTE J1810-197 and PSR J1622-4950.

  3. 3.

Changing $\phi_{0}$. During the revival of radio emission in PSR J1622-4950
(section 4.2 in Camilo et al. 2018), the deflection point of the position
angle is changed horizontally. The overall shape of the position angle changes
from convex to concave. After the magnetar-like activities in the high
magnetic field pulsar PSR J1119-6127, the slope of the position angle even
changes sign (Dai et al. 2018). In our opinion, the twist of the magnetic
field can change the position angle horizontally (equation (27)). This may
explain the change in deflection point in PSR J1622-4950. The problem is that
the changes in the horizontal direction for a twisted dipole field is
relatively small (see discussions following equation (27)). If the field
geometry is more complicated than the twisted dipole case, large changes in
the horizontal direction may be possible.

Furthermore, if there are large changes in the position angle horizontally,
the overall shape of the position angle may change from convex to concave as
can be seen in figure 2. This may correspond to the case of PSR J1622-4950
(Camilo et al. 2018). If the horizontal change of position angle is large
enough, the slope of the position angle can even change sign. This may explain
the observations of PSR J1119-6127 (Dai et al. 2018). For large changes in the
horizontal direction, a more complicated field geometry than the twisted
dipole field is required.

  4. 4.

Wide pulse profile. When the rotating vector model is proposed (Radhakrishnan
& Cooke 1969), it provides an explanation for both the pulse profile and
position angle. The five radio emitting magnetars all have a period about
several seconds. However, their radio pulse profile are very wide, for a duty
cycle about $50\%$. This is much wider than that of normal pulsars (Lyne &
Manchester 1988; Tan et al. 2018). One reason may be that the radio emission
of magnetars originates from a much higher radius in the magnetosphere.
Considering that the magnetars’ magnetic field is twisted, it is more probable
that the polar cap will be larger for a twisted dipole field (Tong 2019).
During untwisting process, the polar cap will be smaller. This may explain the
shrinking hot spot of magnetar X-ray observations (Tong 2019). This may also
explain the wide pulse profile of magnetar radio emissions. One prediction of
the twisted dipole field explanation is that: during the untwisting process,
the pulse profile will become narrower. This is indeed observed during the
diminishing process PSR J1622-4950 (Scholz et al. 2017).

In summary, the ordered position angle of magnetar radio emission tells us
that the large scale magnetic field in the radio emitting region is ordered,
e.g., dipole or twisted dipole. The time varying position angle (change
$\psi_{0}$, slope, $\phi_{0}$) indicates time evolution of the magnetosphere.
A twisted dipole field and its untwisting process may explain the polarization
of magnetar radio emissions. A wide pulse profile and its evolution are also
consistent with a twisted dipole field. More complex field geometry is needed
in order to explain the changes in $\phi_{0}$ quantitatively.

Figure 4: The evolution of $\psi_{0}$ of the radio emitting magnetar PSR
J1622-4950. The data points are from table 1 of Levin et al. (2012) (one data
point with position angle jump is not considered). The solid line is a typical
curve of exponential decay form: $\psi_{0}(t)=-30+25\exp(-(t-t_{0})/2\tau)$
(in unit of degrees), where $t_{0}$ is the date of the first data points, and
the exponential decay timescale is $\tau=90$ days.

##  6 Discussion and conclusion

In order to have high quality polarization measurement, a high flux of the
source is preferred. This means that the polarization measurement of magnetar
is more easily obtained at the beginning of the radio observations (Eatough et
al. 2013). This may also explain why only XTE J1810-197 and PSR J1622-4950
have detailed polarization information (Camilo et al. 2007; Kramer et al.
2007; Dai et al. 2019; Levin et al. 2012; Camilo et al. 2018). The other three
radio emitting magnetars also have polarization observations, but the details
are limited compared with XTE J1810-197 and PSR J1622-4950. Therefore, more
observations are required in order to infer the magnetic field geometry from
radio polarization data.

Considering that there is a continuous transition from normal pulsars to
magnetars, twist of magnetic field may also happen in normal pulsars. This
possibility is strengthened by the two high magnetic field pulsars with
magnetar activities (Gavriil et al. 2008; Archibald et al. 2016). Weak and
strong twist may account for different behaviors of high magnetic field
pulsars and magnetars (Tong 2019). The twist in the case of normal pulsars may
be very weak. And there exist several depolarization effects (Levin et al.
2012). Therefore, the effect of a twisted magnetic field in normal pulsar may
be very difficult to verify. However, with the FAST or SKA telescope, we may
find hints of a twisted magnetic field in some normal pulsars.

In principle, we should fit the observed position angle using the exact
formula (equation (2)): the pure dipole case given by equation (1), the
modification in the twisted dipole case given by equation (21), the twist
parameter $\lambda$ can be either positive or negative, and $\theta_{\rm obs}$
is also determined by the geometry parameters (equation (31)). The twisted
dipole scenario only introduce one additional free parameter $\lambda$ (or
$n$, or $\Delta\phi_{\rm max}$ which are equivallent). From previous examples
of PSR J1622-4950, if a changing $\psi_{0}$ or $\phi_{0}$ or changing slope
etc are seen in one radio emitting magnetar, it may due to an untwisting
magnetosphere. Irrespective of this point, the observed position angle can
also be fitted by the dipole formula. Then the time evolving $\psi_{0}$ can be
modeled by the twisted dipole field etc. This is the second way of fitting the
observations. The merit of the second method is that the changing $\psi_{0}$
can be modeled by the twisted dipole field or other magnetospheric models of
magnetars.

When modeling the X-ray outburst of magnetars, the modeling is often based on
a twisted magnetic field (Thompson et al. 2002; Beloborodov 2009; Pavan et al.
2009; Fujisawa & Kisaka 2014; Akgun et al. 2016; Kojima 2017; Tong 2019). We
suggest that the researchers also calculate the corresponding polarization
position angle in their magnetic field geometry. The procedure can be similar
to that in Section 2 (simple geometry calculation) or Section 3 (more detailed
calculations based on a perturbation analysis). For any magnetic field
geometry, an order of magnitude estimation can be done. According to the
definition of the magnetic field lines (Wolfson 1995)

$\frac{r\sin\theta d\phi}{rd\theta}=\frac{B_{\phi}}{B_{\theta}}.$ (28)

Therefore, in order of magnitude, the changes of the field line in the
$\phi$-direction is

$\Delta\phi_{\rm twist}\sim\frac{B_{\phi}}{\sin\theta
B_{\theta}}\sim\frac{B_{\phi}}{\sin\theta B_{\rm p}},$ (29)

where $B_{\phi}$ is the toroidal magnetic field, and $B_{\rm p}$ is the
poloidal magnetic field. Adding the correct sign, the change in the position
angle is

$\Delta\psi_{\rm twist}\sim-4\frac{B_{\phi}}{\sin\theta B_{\rm p}}.$ (30)

The numerical factor $4$ is obtained by comparison with the detailed
calculations in the twisted dipole case (equation (21)). Therefore, for any
magnetospheric modeling of magnetars, the corresponding modifications to the
position angle can be estimated using equation (30). Only the relative
strength between the toroidal field and the poloidal field is required.

If the toroidal field is due to the polar cap current (Hibschman & Arons
2001), according to equation (30), an order of magnitude estimation is:
$\Delta\psi_{\rm twist}\sim 4(r/R_{\rm L})(J/J_{\rm GJ})\cos\alpha$. Comparing
the estimation here and equation (D13) in Hibschman & Arons (2001), it can be
seen that equation (30) can give the correct sign and correct order of
magnitude estimation for small $\theta$.

In conclusion, the twisted magnetic field of magnetars will bring changes to
their polarization position angle. The relation between the emission point and
the line of sight is slightly modified. Current radio observations of
magnetars may indicate that their magnetic field is indeed twisted. At
present, only the twisted dipole case is calculated. Calculations for more
complicated field geometry (e.g., multipole field) are needed in the future.

## Acknowledgements

H.Tong is supported the Ministry of Science and Technology of the People’s
Republic of China (2020SKA0120300) and NSFC (11773008). The authors would like
to thank R.X.Xu for helpful comments.

## Data availability

This is a theoretical paper, mainly analytical. All the formulae are available
in the article.

## References

  * (1) Akgün T., Miralles J. A., Pons J. A., et al., 2016, MNRAS, 462, 1894
  * (2) Anderson G. E., Gaensler B. M., Slane P. O., et al., 2012, ApJ, 751, 53
  * (3) Archibald R. F., Kaspi V. M., Tendulkar S. P., et al., 2016, ApJ, 829, L21
  * (4) Beloborodov A. M., 2009, ApJ, 703, 1044
  * (5) Blaskiewicz M., Cordes J. M., Wasserman I., 1991, ApJ, 370, 643
  * (6) Camilo F., Ransom S. M., Halpern J. P., et al., 2006, Nature, 442, 892
  * (7) Camilo F., Ransom S. M., Halpern J. P., et al., 2007a, ApJ, 666, L93
  * (8) Camilo F., Reynolds J., Johnston S., et al., 2007b, ApJ, 659, L37
  * (9) Camilo F., Reynolds J., Johnson S., et al., 2008, ApJ, 679, 681
  * (10) Camilo F., Ransom S. M., Halpern J. P., et al., 2016, ApJ, 820, 110
  * (11) Camilo F., Scholz P., Serylak M., et al., 2018, ApJ, 856, 180
  * (12) Carroll B. W., Ostlie D. A., 2018, An introduction to modern astrophysics, Cambridge, Cambridge
  * (13) Coti Zelati F., Rea N., Pons J. A., et al., 2018, MNRAS, 474, 961
  * (14) Dai S., Johnson S., Weltevrede P., et al., 2018, MNRAS, 480, 3584
  * (15) Dai S., Lower M. E., Bailes M., et al., 2019, ApJL, 874, L14
  * (16) Duncan R. C., Thompson C., 1992, ApJL, 392, L9
  * (17) Eatough R. P., Falcke H., Karuppusamy R., et al., 2013, Nature, 501, 391
  * (18) Everett J. E., Weisberg J. M., 2001, ApJ, 553, 341
  * (19) Fujisawa K., Kisaka S., 2014, MNRAS, 445, 2777
  * (20) Gavriil F. P., Gonzalez M. E., Gotthelf E. V., et al., 2008, Science, 319, 1802
  * (21) Griffiths D. J., 2017, Introduction to electrodynamics, Cambridge, Cambridge
  * (22) Hibschman J. A., Arons J., 2001, ApJ, 546, 382
  * (23) Johnston S., Weisberg J. M., 2006, MNRAS, 368, 1856
  * (24) Johnston S., Kramer M., 2019, MNRAS, 490, 4565
  * (25) Kaspi V. M., Beloborodov A. M., 2017, ARAA, 55, 261
  * (26) Kojima Y., 2017, MNRAS, 468, 2011
  * (27) Kramer M., Stappers B. W., Jessner A., et al., 2007, MNRAS, 377, 107
  * (28) Levin L., Bailes M., Bates S., et al., 2010, ApJ, 721, L33
  * (29) Levin L., Bailes M., Bates S. D., et al., 2012, MNRAS, 422, 2489
  * (30) Lower M. E., Shannon R. M., Johnson S., Bailes M., 2020a, ApJL, 896, L37
  * (31) Lower M. E., Johnston S., Shannon R. M., Bailes M., Camilo F., 2020b, arXiv:2011.12463
  * (32) Lyne A. G., Manchester R. N., 1988, MNRAS, 234, 477
  * (33) Manchester R. N., Han J. L., Qiao G. J., 1998, MNRAS, 295, 280
  * (34) Pavan L., Turolla R., Zane S., et al., 2009, MNRAS, 395, 753
  * (35) Radhakrishnan V., Cooke D. J., 1969, ApL, 3, 225
  * (36) Rankin J. M., 1993, ApJ, 405, 285
  * (37) Rea N., Esposito P., Turolla R., et al., 2010, Science, 330, 944
  * (38) Scholz P., Camilo F., Sarkissian J., et al., 2017, ApJ, 841, 126
  * (39) Tan C. M., Bassa C. G., Cooper S., et al., 2018, ApJ, 866, 54
  * (40) Thompson C., Lyutikov M., Kulkarni S. R., 2002, ApJ, 574, 332
  * (41) Tiengo A., Esposito P., Mereghetti S., et al., 2013, Nature, 500, 312
  * (42) Tong H., 2019, MNRAS, 489, 3769
  * (43) Tong H., Huang L., 2020, MNRAS, 497, 2680
  * (44) Wolfson R., 1995, ApJ, 443, 810
  * (45) Weng S. S., Gogus E., Guver T., Lin L., 2015, ApJ, 805, 81
  * (46) Zhou P., Chen Y., Li X. D., et al., 2014, ApJ, 781, L16

##  Appendix A Rotating vector model for a dipole field

The expression of the position angle for a rotating dipole can be found using
spherical trigonometry. From figure 1, the spherical triangle $\triangle RMP$
has two sides and one angle known: $\overline{RM}=\alpha$,
$\overline{RP}=\zeta$, and $\angle MRP=\phi$. The third side
$\overline{PM}=\theta_{\rm obs}$ is the colatitude of the line of sight in the
magnetic frame. It can be found using spherical trigonometry (Carroll & Ostlie
2018). Using the law of cosines for sides

$\cos\theta_{\rm obs}=\cos\alpha\cos\zeta+\sin\alpha\sin\zeta\cos\phi.$ (31)

Using the law of sines, the position angle $\angle MPR=\psi$ is:

$\frac{\sin\alpha}{\sin\psi}=\frac{\sin\theta_{\rm obs}}{\sin\phi}.$ (32)

Therefore, the sine of the position angle is

$\sin\psi=\sin\alpha\frac{\sin\phi}{\sin\theta_{\rm obs}}.$ (33)

Here $\sin\theta_{\rm obs}$ is required. However, only $\cos\theta_{\rm obs}$
is obtained above. Using the law of cosines for sides again for $\alpha$

$\cos\alpha=\cos\zeta\cos\theta_{\rm obs}+\sin\zeta\sin\theta_{\rm
obs}\cos\psi.$ (34)

Then the cosine of the position angle is

$\cos\psi=\frac{\cos\alpha-\cos\zeta\cos\theta_{\rm
obs}}{\sin\zeta\sin\theta_{\rm obs}}.$ (35)

Then it is straight forward to show that

$\tan\psi\equiv\frac{\sin\psi}{\cos\psi}=\frac{\sin\alpha\sin\phi}{\cos\alpha\sin\zeta-\sin\alpha\cos\zeta\cos\phi}.$
(36)

This is the commonly cited formula for the position angle. In practice, $\psi$
should be replaced by $\psi-\psi_{0}$, $\phi$ should be replaced by
$\phi-\phi_{0}$, where $\psi_{0}$ and $\phi_{0}$ are constants to be
determined by fitting the observations (Everett & Weisberg 2001; Johnston &
Kramer 2019).

The third angle $\angle RMP=\phi_{\rm obs}$ is the azimuthal angle of the line
of sight in the magnetic frame. Similar calculations showed that

$\displaystyle\sin\phi_{\rm obs}=\sin\zeta\frac{\sin\phi}{\sin\theta_{\rm
obs}}$ (37) $\displaystyle\cos\phi_{\rm
obs}=\frac{\cos\zeta-\cos\alpha\cos\theta_{\rm obs}}{\sin\alpha\sin\theta_{\rm
obs}}$ (38)

These two equations will be used during the following calculations.

##  Appendix B Detailed calculations for a dipole field

According to Hibschman & Arons (2001), given the magnetic field, the tangent
vector along the magnetic field can be found by normalization
$\bf\hat{t}=B/|B|$. The curvature vector at that point is
${\bf\kappa=(\hat{t}\cdot\nabla)\hat{t}}=1/2\bf\nabla\hat{t}^{2}-\hat{t}\times(\nabla\times\hat{t})$.
The normal vector is the normalized curveture vector
$\bf\hat{n}=\kappa/|\kappa|$. The binormal vector is perpendicular to both the
tangent vector and normal vector $\bf\hat{b}\equiv\hat{t}\times\hat{n}$. The
projection of the rotational axis onto the plane of sky is: $\bf\Omega_{\rm
p}\equiv\hat{\Omega}-(\hat{\Omega}\cdot\hat{t})\hat{t}$, where
$\bf\hat{\Omega}$ is the unit vector along the roational axis. The angle
between the rotational axis and the normal vector is the position angle, shown
in figure 5. The angle between the rotation axis and the binormal vector is:
$\cos\psi^{\prime}=\bf\hat{b}\cdot\hat{\Omega}_{\rm p}$,
$\sin\psi^{\prime}=\bf\hat{t}\cdot(\hat{b}\times\hat{\Omega}_{\rm p})$.
$\bf\hat{\Omega}_{\rm p}$ is the unit vector along $\bf\Omega_{\rm p}$.
Therefore, the position angle is

$\displaystyle\tan\psi$ $\displaystyle=$
$\displaystyle-\cot\psi^{\prime}=\frac{\bf\hat{b}\cdot\Omega_{\rm
p}}{\bf\hat{b}\cdot(\hat{t}\times\Omega_{\rm p})}$ (39) $\displaystyle=$
$\displaystyle\frac{\bf\hat{b}\cdot\hat{\Omega}}{\bf\hat{b}\cdot(\hat{t}\times\hat{\Omega})},$
(40)

where during the last step further simplification are used.

For a dipole magnetic field ${\bf
B_{0}}=(\frac{2\mu\cos\theta}{r^{3}},\frac{\mu\sin\theta}{r^{3}},0)$, it can
be shown that the final position angle is the same as in equation (1). The
power of this detailed calculation process will be manifested when there are
perturbations to the magnetic field (Hibschman & Arons 2001). The magnetic
field may be rewritten as $\bf B=B_{0}+B_{1}$, where $\bf B_{0}$ is the dipole
magnetic field, $\bf B_{1}$ is the perturbed magnetic field. Then the tangent
vector $\bf\hat{t}=\hat{t}_{0}+t_{1}$, curvature
$\bf\kappa=\kappa_{0}+\kappa_{1}$, normal vector
$\bf\hat{n}=\hat{n}_{0}+n_{1}$, and binormal vector
$\bf\hat{b}=\hat{b}_{0}+b_{1}$ will all be perturbed. Keeping only the linear
term, the perturbed quantities are

$\displaystyle\bf t_{1}$ $\displaystyle=$ $\displaystyle\bf B_{1}/|B_{0}|$
(41) $\displaystyle\bf\kappa_{1}$ $\displaystyle=$
$\displaystyle\bf(\hat{t}_{0}\cdot\nabla)t_{1}+(t_{1}\cdot\nabla)\hat{t}_{0}$
(42) $\displaystyle=$ $\displaystyle\bf\nabla(\hat{t}_{0}\cdot
t_{1})-\hat{t}_{0}\times(\nabla\times
t_{1})-t_{1}\times(\nabla\times\hat{t}_{0})$ (43) $\displaystyle\bf n_{1}$
$\displaystyle=$
$\displaystyle\bf\frac{1}{|\kappa_{0}|}(\kappa_{1}-(\hat{n}_{0}\cdot\kappa_{1})\hat{n}_{0})$
(44) $\displaystyle\bf b_{1}$ $\displaystyle=$
$\displaystyle\bf\hat{t}_{0}\times n_{1}+t_{1}\times\hat{n}_{0}.$ (45)

The modification of position angle is

$\Delta\psi=\frac{\bf\hat{t}_{0}\cdot(b_{1}\times\Omega_{\rm
p})}{\bf\hat{b}_{0}\cdot\Omega_{\rm p}}.$ (46)

In Hibschman & Arons (2001), they also considered the change of binormal
vector due to changes in the emission point. During this process, more
assumptions are adopted (e.g., the radio emission arises at a given radius,
only lateral motion is allowed). Furthermore, different perturbations can be
simply be added since only linear terms are kept. This effect (changes in
emission point) and others will not affect the changes due to a perturbed
magnetic field. The final result will be simple addition of different effects.
Therefore, in the main text we do not consider changes in the emission point
and only consider the change of position angle due to a twisted magnetic
field.

Figure 5: Position angle between the normal vector and the projection of
rotational axis.

##  Appendix C The relation between emission point and line of sight for a
dipole field

The line of sight lies in the direction ($\theta_{\rm obs}$, $\phi_{\rm
obs}$), in spherical coordinate, in the magnetic frame. In Cartesian
coordinate

$\bf\hat{l}=\sin\theta_{\rm obs}\cos\phi_{\rm obs}\hat{x}+\sin\theta_{\rm
obs}\sin\phi_{\rm obs}\hat{y}+\cos\theta_{\rm obs}\hat{z}.$ (47)

For a dipole field, the tangent vector at ($\theta$, $\phi$) is

$\displaystyle{\bf\hat{t}_{0}}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{1}}(\cos\theta,\frac{1}{2}\sin\theta,0)$ (48)
$\displaystyle=$
$\displaystyle\frac{1}{N_{1}}(\cos\theta\hat{r}+\frac{1}{2}\sin\theta\hat{\theta}),$
(49)

where $N_{1}=\sqrt{1-(3/4)\sin^{2}\theta}$ is the normalization factor. The
relation between the base vector in spherical coordinate and Cartesian
coordinate is (Griffiths 2017)

$\displaystyle\hat{r}$ $\displaystyle=$
$\displaystyle\sin\theta\cos\phi\hat{x}+\sin\theta\sin\phi\hat{y}+\cos\theta\hat{z}$
(50) $\displaystyle\hat{\theta}$ $\displaystyle=$
$\displaystyle\cos\theta\cos\phi\hat{x}+\cos\theta\sin\phi\hat{y}-\sin\theta\hat{z}$
(51) $\displaystyle\hat{\phi}$ $\displaystyle=$
$\displaystyle-\sin\phi\hat{x}+\cos\phi\hat{y}.$ (52)

Then, the tangent vector rewritten Cartesian coordinate is

${\bf\hat{t}_{0}}=\frac{1}{N_{1}}[\frac{3}{2}\sin\theta\cos\theta\cos\phi\hat{x}+\frac{3}{2}\sin\theta\cos\theta\sin\phi\hat{y}+(\cos^{2}\theta-\frac{1}{2}\sin^{2}\theta)\hat{z}]$
(53)

It is a general assumption that the tangent vector should be in the direction
of the line of sight (Hibschman & Arons 2001). By equaling the tangent vector
and the line of sight unit vector, it is found that

$\displaystyle\frac{1}{N_{1}}\frac{3}{2}\sin\theta\cos\theta\cos\phi$
$\displaystyle=$ $\displaystyle\sin\theta_{\rm obs}\cos\phi_{\rm obs}$ (54)
$\displaystyle\frac{1}{N_{1}}\frac{3}{2}\sin\theta\cos\theta\sin\phi$
$\displaystyle=$ $\displaystyle\sin\theta_{\rm obs}\sin\phi_{\rm obs}$ (55)
$\displaystyle\frac{1}{N_{1}}(\cos^{2}\theta-\frac{1}{2}\sin^{2}\theta)$
$\displaystyle=$ $\displaystyle\cos\theta_{\rm obs}.$ (56)

From the first two equations, it can be inferred that $\phi=\phi_{\rm obs}$.
Then both the sine and cosine of $\theta_{\rm obs}$ is known. It can combined
to obtain

$\tan\theta_{\rm
obs}=\frac{3}{2}\frac{\tan\theta}{1-\frac{1}{2}\tan^{2}\theta}.$ (57)

This a quadratic equation of $\tan\theta$ and can be solved analytically. For
the case of small $\theta$, it can be found directly that (Hibschman & Arons
2001)

$\theta_{\rm obs}=\frac{3}{2}\theta.$ (58)


2022

These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Vivek Muthurangu

[1]UCL Institute of Cardiovascular Science, University College London, London,
United Kingdom

2]Great Ormond Street Hospital, London, United Kingdom

3]Children’s Hospital at Westmead, Sydney, Australia

4]Faculty of Medicine and Health, University of Sydney, Sydney, Australia

# Automatic Segmentation of the Great Arteries for Computational Hemodynamic
Assessment

Javier Montalt-Tordera<EMAIL_ADDRESS>Endrit Pajaziti
<EMAIL_ADDRESS>Rod Jones<EMAIL_ADDRESS>Emilie Sauvage
<EMAIL_ADDRESS>Rajesh Puranik<EMAIL_ADDRESS>Aakansha Ajay
Vir Singh<EMAIL_ADDRESS>Claudio Capelli
<EMAIL_ADDRESS>Jennifer Steeden<EMAIL_ADDRESS>Silvia
Schievano<EMAIL_ADDRESS><EMAIL_ADDRESS>*  [  [  [

###### Abstract

Background: Computational fluid dynamics (CFD) is increasingly used for the
assessment of blood flow conditions in patients with congenital heart disease
(CHD). This requires patient-specific anatomy, typically obtained from
segmented 3D cardiovascular magnetic resonance (CMR) images. However,
segmentation is time-consuming and requires expert input. This study aims to
develop and validate a machine learning (ML) method for segmentation of the
aorta and pulmonary arteries for CFD studies.

Methods: 90 CHD patients were retrospectively selected for this study. 3D CMR
images were manually segmented to obtain ground-truth (GT) background, aorta
and pulmonary artery labels. These were used to train and optimize a U-Net
model, using a 70-10-10 train-validation-test split. Segmentation performance
was primarily evaluated using Dice score. CFD simulations were set up from GT
and ML segmentations using a semi-automatic meshing and simulation pipeline.
Mean pressure and velocity fields across 99 planes along the vessel
centrelines were extracted, and a mean average percentage error (MAPE) was
calculated for each vessel pair (ML vs GT). A secondary observer (SO)
segmented the test dataset for assessment of inter-observer variability.
Friedman tests were used to compare ML vs GT, SO vs GT and ML vs SO metrics,
and pressure / velocity field errors.

Results: The network’s Dice score (ML vs GT) was 0.945 (interquartile range:
0.929–0.955) for the aorta and 0.885 (0.851–0.899) for the pulmonary arteries.
Differences with the inter-observer Dice score (SO vs GT) and ML vs SO Dice
scores were not statistically significant for either aorta or pulmonary
arteries ($p=0.741$, $p=0.061$). The ML vs GT MAPEs for pressure and velocity
in the aorta were 10.1% (8.5–15.7%) and 4.1% (3.1–6.9%) respectively, and for
the pulmonary arteries 14.6% (11.5–23.2%) and 6.3% (4.3–7.9%), respectively.
Inter-observer (SO vs GT) and ML vs SO pressure and velocity MAPEs were of a
similar magnitude to ML vs GT ($p>0.2$).

Conclusions: ML can successfully segment the great vessels for CFD, with
errors similar to inter-observer variability. This fast, automatic method
reduces the time and effort needed for CFD analysis, making it more attractive
for routine clinical use.

###### keywords:

segmentation, machine learning, neural network, computational fluid dynamics,
congenital heart disease, magnetic resonance imaging

##  1 Background

The past two decades have seen increasing interest in the use of computational
fluid dynamics (CFD) for the assessment of cardiovascular disease, including
congenital heart disease (CHD) Biglino2017ComputationalTranslation .
Computational fluid dynamics models enable realistic calculation of patient-
specific blood flow conditions and provide valuable insights into pathological
hemodynamics. These models can also be used to predict hemodynamic response to
interventions, thereby aiding therapeutic planning.

The patient specific anatomies needed for CFD are often derived from three
dimensional (3D) cardiovascular magnetic resonance (CMR), particularly cardiac
and respiratory gated whole-heart sequences. This is because whole-heart
images have the sharp borders and high contrast necessary for semi-automated
segmentation of cardiovascular structures.

Automatic or semi-automatic segmentation methods based on shape models or
level-set algorithms have existed for years Wang2014FastPropagation ;
Weese2001ShapeSegmentation ; Peters2007AutomaticVolumes . However, they often
require some type of user input or manual post-correction, rely on priors
which are not readily available, or may struggle to adapt to abnormal
anatomies (e.g., CHD). Thus, segmentation remains one of the most user-
intensive and time-consuming parts of the CFD workflow, and one of the
barriers to greater clinical use.

Recently, it has been shown that machine learning (ML) can accurately segment
ventricles and great vessels from CMR images Chen2020DeepReview ;
Berhane2020FullyLearning ; Zhuang2019EvaluationChallenge . Quantitative
metrics derived from ML segmentations (e.g., ventricular volumes) compare well
with manual segmentations Bai2018AutomatedNetworks ;
Bottcher2020FullyAlgorithm  and these techniques are now entering clinical
practice. However, the effectiveness of ML segmentation for CFD has not
previously been investigated.

The aims of this study were to: (i) develop a ML method for simultaneous
segmentation of the aorta and pulmonary arteries from whole heart CMR images
in patients with paediatric or adult CHD, (ii) compare conventional and ML
segmentations using traditional image-based scores, (iii) compare CFD metrics
derived from both conventional and ML segmentations, and (iv) investigate the
association between image-based scores and CFD errors.

##  2 Methods

###  2.1 Subjects

90 cardiac triggered, respiratory navigated, 3D whole heart, balanced, steady
state, free precession (WH-bSSFP) data were collected from previously scanned
children and adults with paediatric or congenital heart disease (excluding
patients with single ventricles). All patients were scanned on a 1.5T scanner
(Avanto, Siemens Healthineers AG, Erlangen, Germany) using a standard WH-bSSFP
sequence Steeden2020RapidSuper-resolution . The imaging protocol was as
follows: orientation: sagittal, matrix size: 256 x 144 x 96 (head-foot,
anterior-posterior, left-right), acquired voxel size: 1.6 mm (isotropic), flip
angle: 90°. Image acquisition was accelerated using GRAPPA (factor of 2 along
phase encoding dimension) and partial Fourier (factor of 6/8 along both phase
and slice encoding dimensions). The use of retrospectively collected training
and test data was approved by the local research ethics committee, and written
consent was obtained from all subjects/guardians (Ref: 06/Q0508/124).

Additionally, 10 external examples were retrospectively collected from a
different centre. These were scanned on a 1.5 T scanner (Ingenia, Philips
Healthcare, Amsterdam, Netherlands) with the following imaging protocol:
orientation: axial, matrix size: 240 x 240 x 110 (left-right, anterior-
posterior, head-foot), acquired voxel size: 1.44 mm (isotropic), flip angle:
90°. Image acquisition was accelerated using SENSE (reduction factor of 2) and
partial Fourier (6/8). Collection of this data and sharing with our site was
approved by the local research ethics committee, and written consent was
obtained from all subjects/guardians (Protocol No. X20-0237 and
2020/ETH01333).

###  2.2 Ground Truth Segmentation

Reference standard conventional segmentation of the aorta and pulmonary
arteries was performed using a semi-automatic technique with manual correction
(Plug-ins built in Horos v4.0, Horosproject.org sponsored by Nimble Co LLC
d/b/a Purview, Maryland, USA). Initial segmentation was done using the fast
level-set method Wang2014FastPropagation . This requires the user to: (i) set
a threshold, (ii) place seeds in the vessel of interest and (iii) add blocking
regions to prevent segmentation of unwanted structures. The quality of this
initial segmentation is dependent on both the underlying anatomy and the image
quality, but manual correction is always required to remove unwanted
structures and clip vessels. The proximal limit of both the aortic and
pulmonary artery segmentations was the semi-lunar valve. The distal limit of
the segmentations were the diaphragmatic level of the aorta, and hilar
branches of the pulmonary arteries. Head and neck arteries were manually
removed at their origin.

All 90 datasets were segmented by a primary observer (RJ — 10 years’
experience in CMR post-processing). We refer to the primary observer’s
segmentations as the ground truth (GT). In addition, a secondary observer (VM
— 19 years’ experience in CMR post-processing) segmented 10 of these images
(test set, see below) to investigate inter-observer variability. We refer to
these as the secondary observer (SO) data. The 10 external examples were
segmented by the secondary observer only following the same procedure.

###  2.3 Data Preparation

Prior to ML training, the pixel intensities of the WH-bSSFP data were
normalized (range [0, 1]). The aortic and pulmonary binary segmentation masks
were concatenated in the channel dimension and combined with a third channel
containing a binary “background” mask (one-hot encoding). The images and
3-channel segmentation masks were either centrally cropped or symmetrically
zero-padded to a fixed matrix size of 160×96×64 (superior-inferior, anterior-
posterior, left-right). Finally, the image-label pairs were randomly split
into a training set (70 examples, 78%), a validation set (10 examples, 11%)
and a test set (10 examples, 11%). This split was used to maximize the size of
the training set, while providing sufficient data for validation and
statistical analysis.

The examples from the external test set were reoriented, interpolated and
cropped to match the orientation, matrix size and voxel size of the in-house
data.

###  2.4 Network Architecture

A U-Net Ronneberger2015U-Net:Segmentation  convolutional neural network was
used to simultaneously segment the aorta and pulmonary arteries from WH-bSSFP
data. The network architecture is shown in Fig. 1. Each convolutional layer
was followed by a batch normalization layer and a rectified linear unit (ReLU)
activation. Downscaling was performed using max-pooling layers and upscaling
was performed using transpose convolution layers. The number of convolutional
filters after the first layer was set to double after each downscaling layer
and halve after each upscaling step. The final convolutional layer has three
filters (equalling the number of possible classes – aorta, pulmonary artery
and background), followed by a softmax activation. Final predicted labels were
obtained by assigning each pixel to the class with the highest probability.

Figure 1: Parameterized network architecture. The height of the blocks
represents changes in spatial resolution while the width represents the number
of filters or channels. The figure shows the three network parameters whose
value was optimized: initial filters, layers per block and scales. This
example is shown with 2 layers per block and 3 scales. BN: batch
normalization.

###  2.5 Training, Hyperparameter Optimization and Evaluation

The network implementation and training scheme were parametrized to allow
investigation of multiple hyperparameter values, with the full search space
shown in Table 1. We used the Hyperband algorithm Li2018Hyperband:Optimization
to perform efficient hyperparameter optimization. This method samples the
search space randomly and adaptively allocates more computational resources to
the most promising hyperparameters combinations. Dice score was used by the
Hyperband to assess performance and choose the final model.

Table 1: Hyperparameter search space.

Parameter | Type | Domain
---|---|---
Scales | Architecture | {2, 3, 4}
Layers per block | Architecture | {2, 3, 4}
Initial filters | Architecture | {32, 64}
Learning rate | Training | [0.0001, 0.01]111The learning rate was sampled using a log-uniform distribution.
Batch size | Training | {2, 4}
Loss function | Training | {CCE222CCE: categorical cross-entropy., Dice, IoU, Tversky, focal Tversky}
|  |

00footnotetext: $\\{\cdot\\}$ represents a discrete set of values; $[\cdot]$
represents a continuous interval of real values.

The neural network and related functionality were implemented and trained
using TensorFlow Abadi2016TensorFlow . In particular, the model
implementation, losses and metrics are available in TensorFlow MRI
TensorFlowMRI , an open-source framework developed in-house. Weights were
initialized using He’s method He2015DelvingClassification  and optimized using
the Adam algorithm Kingma2015Adam:Optimization . Training, including
hyperparameter optimization, took 24 hr on an Nvidia Titan RTX GPU with 24 GB
of onboard RAM (Nvidia Corporation, Santa Clara, CA, USA).

The optimized ML model was evaluated on the test dataset against the ground
truth segmentations (ML vs GT). The accuracy of segmentation was quantified
using several image-based segmentation metrics: Dice score, Intersection over
Union (IoU), Hausdorff distance (HD) and average surface distance (ASD). Each
metric was computed independently for each vessel. Additionally, the same
metrics were calculated for the secondary observer’s segmentation against the
GT (SO vs GT), and between the ML model and the secondary observer (ML vs SO).

To assess generalization ability to data from other sources, the ML model was
also evaluated on the external test set. We computed the same set of metrics
(Dice, IoU, HD and ASD) between the ML predictions (Ext-ML) and the ground
truth segmentations (we refer to these as Ext-SO, because the manual
segmentation was performed by the same person as the SO data).

Prior to evaluation, ML masks were filtered to remove all but the largest
connected component, as identified using 3D connected component labelling with
26-connectivity Silversmith2021Cc3d:Images . This postprocessing step was used
to eliminate small background regions which had been misclassified as vessels.

###  2.6 Surface and Volume Meshing Pipeline

The resultant segmentation masks were converted into finite element volume
meshes, using the processes shown in Fig. 2. The masks (GT, ML and SO) were
first transformed into surface meshes, by applying the marching cubes
algorithm (implemented in VMTK), and then re-meshed and smoothed with
consistent parameters. The surface meshes were clipped manually at inlets and
outlets to create planar surfaces; 30 mm extensions were added to ensure the
fluid flow in the region of interest was fully developed and capped at the
ends to generate close surfaces. These were meshed with tetrahedral elements
to build the final unstructured grid for CFD analysis. The grid resolution was
determined through a sensitivity analysis (see Appendix A).

To assess the effect of the manual clipping of the anatomies on the CFD, the
simulations from the ML segmentations were also re-run using the same inlet
and outlet locations as the GT data, defined by overlaying the GT clipped
geometry to the corresponding ML geometry.

Figure 2: Automatic mesh processing pipeline from segmentation to CFD
analysis, followed by post-processing to reshape the data in a consistent
format between subjects (99 planes from inlet to outlet containing average
pressure and velocity).

###  2.7 Computational Fluid Dynamics and Boundary Conditions

CFD simulations were carried out using the solver Fluent 19.0 (Ansys,
Canonsburg, PA, USA). Blood was modelled as an incompressible Newtonian fluid
with density 1060 kg/m3 and 0.004 Pa$\cdot$s dynamic viscosity
Gijsen1999TheModel . Vessel walls were considered rigid, and no-slip
conditions were imposed. A laminar, steady-state model was selected to
simulate blood flow at peak systole Caballero2013AAorta ;
Bonfanti2017ComputationalData . A generalisable inlet condition for the aorta
and pulmonary artery was applied to all subjects with a uniform (plug) inlet
velocity profile of 0.66 m/s for the aorta and 0.57 m/s for the pulmonary
artery Gabe1969MeasurementProbe . The outlets for all cases were assumed to be
at zero pressure and the convergence criteria was set at $10^{-4}$ for the
residual errors. Simulations were run on a Dell workstation, with a Xeon CPU
E5-2630 (24 processors at 2.3 GHz), 32 GB RAM and an Nvidia GeForce GTX 1080
Ti.

###  2.8 CFD Post-Processing and Analysis

To compare the flow field between the pairs of different unstructured meshes
(ML vs GT, SO vs GT and ML vs SO), correspondence was created by subdividing
each vessel with 99 planes orthogonal to the centrelines, calculated in VMTK,
and equally distanced. Static pressure and velocity magnitude were averaged in
each plane (see Fig. 2) and a percentage error was calculated for each plane
pair, using ML as reference (GT in the comparison between SO and GT). The mean
absolute percentage errors (MAPE) for pressure and velocity were computed for
each vessel pair.

###  2.9 Statistical Analysis

Shapiro-Wilk tests were used to test the normality of the different
segmentation metrics and CFD errors, grouped by vessel (aorta and pulmonary
artery), and segmentation pair (ML vs GT, ML vs SO and SO vs GT). Wilcoxon
signed rank tests were used to compare the pressure and velocity errors for
the ML vs GT group. Mann-Whitney U-tests were used to compare segmentation
metrics and flow field errors between the aorta and the pulmonary artery, for
the ML vs GT group. Friedman tests for repeated measurements were performed to
compare segmentation metrics and flow field errors between the ML vs GT, ML vs
SO and SO vs GT groups, for both aorta and pulmonary artery segmentations.
Significant Friedman test results were followed up by pairwise Wilcoxon post-
hoc tests. Additionally, Wilcoxon signed rank tests were used to compare ML vs
GT and SO vs GT metrics for both aorta and pulmonary artery segmentations.
Mann-Whitney U-tests were used to compare Ext-ML vs Ext-SO segmentation
metrics against ML vs SO metrics. Wilcoxon signed rank tests were used to
compare the pressure and velocity errors for the manually clipped ML vs GT
data and the equally clipped ML vs GT data. Pearson’s correlation coefficient
was used to measure the linear relationship between each pair of a
segmentation metric (i.e., Dice, IoU, HD or ASD) and a flow field error
(pressure or velocity MAPEs), for both aorta and pulmonary artery
segmentations. The _p_ -value was calculated for each comparison to test non-
correlation. Throughout this work, a _p_ -value $<$ 0.05 was considered
statistically significant.

##  3 Results

###  3.1 Hyperparameter Optimization

A total of 124 hyperparameter configurations were sampled during the neural
network optimization procedure (see Appendix B). The best performing
configuration was as follows: scales = 3, layers per block = 2, initial
filters = 64, learning rate = 3.46x10-4, batch size = 2, and loss function =
focal Tversky. This model was selected and used in all further experiments.

###  3.2 ML Segmentation

The ML segmentation was successful in all 10 test datasets. The specific
diagnoses for these patients were: repaired tetralogy of Fallot ($n=1$),
repaired Tetralogy of Fallot with mild right pulmonary artery stenosis
($n=1$), Marfan syndrome with dilated aorta ($n=1$), Marfan syndrome with
pectus excavatum ($n=1$), dilated pulmonary artery ($n=1$), bicuspid aortic
valve with dilated aorta and unrepaired VSD ($n=1$), repaired double outlet
right ventricle with right sided arch ($n=1$), unrepaired atrial septal defect
($n=1$), aortic regurgitation with dilated aorta ($n=1$), post Ross procedure
with mechanical aortic valve ($n=1$). Inference time for the ML model was
approximately 160 ms for simultaneous segmentation of aorta and pulmonary
arteries (compared to approximately 30 minutes for manual segmentation of
aorta and pulmonary arteries). There was good agreement between the ML and GT
segmentation with a median Dice score of 0.945 (interquartile range:
0.929–0.955) for the aorta and 0.885 (0.851–0.899) for the pulmonary arteries.
The Dice score was significantly higher for the aorta than the pulmonary
arteries ($p=0.002$) with similar findings observed for IoU, HD and ASD (Fig.
3A–D).

Figure 3: Segmentation metrics and flow field errors. Three segmentations are
compared in a pairwise fashion: machine learning (ML), ground truth (GT) and
secondary observer (SO). (A–B): Confusion-based similarity metrics: Dice score
and IoU. (C–D) Distance-based similarity metrics: Hausdorff distance (HD) and
average surface distance (ASD), measured in pixels. (E–F) CFD-derived pressure
and velocity mean average percentage errors (MAPE).

The best, median and worst segmented images in terms of Dice score are shown
in Fig. 4. The three main differences were: (i) the length of the vessel
segmented, (ii) differences in pixel labelling that resulted in small
deviations of the vessel border, and (iii) small protrusions at origin of the
carotid and subclavian arteries in the ML segmentations of the aorta.

Figure 4: Test set segmentation overlays. Predicted and ground truth masks are
overlayed over the original images for the best, median and worst test cases.
Aorta and pulmonary artery masks are shown in red and blue, respectively.
Multiplanar reformats of the original 3D volume were manually selected on a
case-by-case basis to be most informative. Best case had Ross procedure and
mechanical aortic valve, the median case had an atrial septal defect and the
worst case had a dilated pulmonary artery.

The aortic inter-observer Dice score (SO vs GT) was 0.949 (0.916–0.960) and
was not significantly different from ML vs GT ($p=0.575$). The pulmonary Dice
score for the SO vs GT was 0.882 (0.870–0.894) and was also not significantly
different from ML vs GT ($p=0.721$). The ML vs SO Dice score was 0.933
(0.924–0.944) for the aorta, which was not significantly different from ML vs
GT and SO vs GT ($p=0.741$), and 0.843 (0.791–0.860) for the pulmonary
arteries, which trended towards being lower than ML vs GT and SO vs GT
($p=0.061$).

The ML segmentation was also successful in the external dataset. The specific
diagnoses for these patients were: cardiomyopathy ($n=4$), normal anatomy
($n=1$), repaired tetralogy of Fallot ($n=1$), left pulmonary artery stenosis
($n=1$), anomalous pulmonary venous drainage ($n=1$), repaired coarctation of
the aorta with hypoplastic arch ($n=1$), bicuspid aortic valve with severe AR
and dilated aortic root ($n=1$). The best, median and worst examples from the
external test set are shown in Fig. 5. There was reasonable agreement between
the Ext-ML and Ext-SO segmentations, with a median Dice score of 0.913
(0.889–0.927) for the aorta and 0.751 (0.728–0.797) for the pulmonary
arteries. Agreement was significantly lower than ML vs SO for the pulmonary
arteries ($p=0.011$), but not for the aorta ($p=0.089$). Similar findings were
observed for IoU, HD and ASD (Fig. 6).

Figure 5: External test set segmentation overlays. Predicted and ground truth
masks are overlayed over the original images for the best, median and worst
test cases. Aorta and pulmonary artery masks are shown in red and blue,
respectively. Multiplanar reformats of the original 3D volume were manually
selected on a case-by-case basis to be most informative. Best case had normal
anatomy with gothic arch, the median case had bicuspid aortic valve and the
worst case had mild left pulmonary artery stenosis. Figure 6: Segmentation
metrics. The model’s segmentation (ML) is compared with the an observer’s
segmentation (SO) for two different datasets: our original data (ML vs SO) and
an external test set from a different site and vendor (Ext-ML vs Ext-SO).
(A–B): Confusion-based similarity metrics: Dice score and intersection-over-
union (IoU). (C–D) Distance-based similarity metrics: Hausdorff distance (HD)
and average surface distance (ASD), measured in pixels.

###  3.3 CFD Metrics

There was overall good agreement in CFD metrics calculated using ML and GT
segmentations (Fig. 3E–F). The median MAPE for pressure and velocity in the
aorta were 10.1% (interquartile range: 8.5–15.7%) and 4.1% (3.1–6.9%)
respectively, and for the pulmonary arteries 14.6% (11.5–23.2%) and 6.3%
(4.3–7.9%). Pulmonary artery MAPEs trended towards higher values compared to
aortic MAPEs, but this did not reach statistical significance ($p=0.081$ for
pressure and $p=0.093$ for velocity). However, pressure was more sensitive
than velocity to different segmentations, with pressure MAPE being 2.5x
greater than velocity MAPE ($p<0.001$).

Fig. 7 shows the surface meshes of test cases with the highest and lowest CFD
MAPE, as well as pressure and velocities along the length of each vessel. Fig.
8 shows pressure and velocity fields calculated using both ML and GT manual
segmentations. The main difference in the surface meshes (particularly for the
worst cases) were associated with the inlets and outlets (angle and size) and
these differences propagated into pressure and velocity fields.

Figure 7: Best and worst cases for aortas and PAs. Graphs show planar-averaged
metrics along the length of the vessels, starting from the inlet. Black
geometries correspond to ground truth, whereas red correspond to predictions.
Best aorta case had atrial septal defect and the worst had a repaired double
outlet right ventricle with right arch. The best PA has a bicuspid aortic
valve with dilated aorta and unrepaired VSD and the worst is a Marfan syndrome
with dilated aorta. Figure 8: Best and worst aorta and pulmonary artery
predictions. Flow fields of pressure and velocity displayed.

SO vs GT (inter-observer) and ML vs SO pressure and velocity MAPEs were of a
similar magnitude to the errors from the ML segmentations (Fig. 3E–F,
$p>0.2$). When the clipping planes of the GT segmentations were used on the ML
geometries, the median pressure and velocity MAPEs were reduced to 8.0/3.1%
($p<0.01$) for the aorta, and to 10.4/3.7% for the pulmonary artery ($p<0.01$)
(see Appendix C).

Fig. 9 illustrates the relationship between the segmentation metrics and the
CFD errors on the ML vs GT comparison. No statistically significant
correlations were found between any of the metrics, either for the aorta or
the pulmonary arteries, for either manual or equally clipped data.

Figure 9: Flow errors against similarity metrics. The figure shows a scatter-
plot matrix where each point corresponds to a subject. In the abscissas, two
confusion-based metrics, Dice and IoU, and two distance-based metrics, the
Hausdorff distance and the average surface distance, measured in pixels. In
the ordinates, the pressure and velocity mean average percentage errors
(MAPE). All values are for the ML vs GT comparison. Red and blue colours
identify aorta and pulmonary artery data, respectively. Trend lines are least-
squares polynomial fits of degree 1. For Dice and IoU, higher is better (more
similar). For Hausdorff distance, average surface distance and pressure and
velocity MAPEs, lower is better.

##  4 Discussion

In this study, a deep neural network was trained to simultaneously segment the
aorta and pulmonary arteries from 3D MRI data. As its primary purpose was to
provide patient specific anatomies for CFD models, we evaluated accuracy using
conventional image-based segmentation metrics and resulting errors in CFD
measures. The main findings were: (i) the proposed network achieved high
performance in terms of image-based segmentation metrics, (ii) there was
reasonable agreement between CFD models derived from the ML and ground truth
manual segmentation, (iii) these errors were similar in magnitude to those
observed between two different manual segmentations, and (iv) there was no
relationship between the segmentation metrics and the resulting CFD errors.

###  4.1 ML Segmentation

In data from the same distribution as the training data, the segmentation
model achieved comparable or better performance than previously reported 3D ML
segmentation of the great vessels, including in patients with congenital heart
disease Berhane2020FullyLearning ; Payer2017Multi-labelConfigurations . This
suggests that the chosen network architecture and subsequent hyperparameter
optimization were sufficient for accurate segmentation. Nevertheless, there
were some differences between the GT and ML segmentations and visual
inspection reveals three main types of error. The first error was a tendency
for ML to start and stop segmenting at slightly different points in the vessel
compared to the ground truth. The second type of error was the presence of
“bumps”, due to the segmentation masks bleeding out at the locations of
arterial branches, particularly in the aorta. Both these errors can be
considered failures to properly demarcate vessel limits, rather than failures
to correctly label blood pool pixels. The third type of error was inaccurate
labelling of blood pixels at the vessel border, resulting in subtle
differences in surface geometry. It should be noted that none of these
patients had very abnormal pulmonary vascular or aortic anatomy, which was
necessary to ensure that CFD models could be created from the segmentations.
However, further testing on complex congenital heart disease is necessary if
segmentation models are to be used more widely. Extension to complex disease
may require further enhancements, and several strategies could potentially
help improve the ML segmentation accuracy and generalizability. These include
increasing the amount and heterogeneity of training data, or performing data
augmentation, both of which improve generalizability and performance of ML
models Sun2017RevisitingEra ; Krizhevsky2012ImageNetNetworks . Another
interesting option might be the inclusion of statistical shape models
Bhalodia2018DeepSSM:Images ; Raju2021DeepDelineation , which could help ensure
that the segmented shapes conform to common patterns.

The model was also tested on 3D data acquired on a different vendor scanner.
Although the type of sequence (3D whole-heart bSSFP) and imaging protocol were
similar to the original data, there were visually apparent differences in
image quality and characteristics. Nevertheless, we observed reasonable
segmentation quality. For the aorta, agreement with a human observer was only
slightly lower than agreement with the same observer in our original data.
However, there was a larger reduction in agreement for the pulmonary arteries.
This suggests there is scope for improving the generalizability of the model.
One of the best solutions for this is to include multi-site, multi-vendor data
in the training set, but this would incur obvious labelling costs and
potential data sharing difficulties. Other approaches to improve robustness to
out-of-distribution data might be the use of data augmentation (e.g., domain
translation methods to generate multi-vendor datasets Zhu2017UnpairedNetworks
; Yan2019TheUnet-GAN ) and the use of strategies that incorporate additional
domain knowledge Liu2021Anatomy-aidedReview .

Of course, segmentation is a challenging task and we demonstrated that the
agreement between two humans was similar to the agreement between ML and the
GT human segmentation. This suggests that ML “errors” are approximately at the
level of the inter-observer variability and similar observations have
previously been made for aortic segmentation Berhane2020FullyLearning . Thus,
we believe ML can provide segmentation with ‘real world’ accuracy.
Furthermore, there are significant advantages of ML over manual segmentation
including very fast segmentation without user interaction and perfect
reproducibility, due to its deterministic nature. This makes ML particularly
useful for removing clinical bottlenecks and accelerating population-based
research.

###  4.2 Relationship Between CFD and Segmentation Errors

We demonstrated reasonable agreement in velocity and pressure fields
calculated from ML and manual segmentations. Importantly, the differences in
CFD metrics using ML vs manual segmentations were of a similar magnitude to
those between two independent manual segmentations. This suggests that ML can
be successfully used to provide a starting point for CFD simulations, with
accuracy similar to inter-observer variability. However, there were some
differences in CFD metrics between ML vs GT segmentations, particularly for
pressure calculations. We think pressure errors are higher because local
deviations in surface geometry tend to cause only local velocity field
derangement, but have a global effect on upstream pressures. This can be seen
in the worst-case aorta, where a kink in the ground truth descending aorta
results in localized flow acceleration, and significantly altered upstream
pressures.

Interestingly, we found no significant correlations between image-based
segmentation metrics and errors in the pressure and velocity fields. This
suggests that neither overlap-based (Dice, IoU) nor boundary distance-based
(HD, ASD) metrics can accurately capture the features that ensure CFD
accuracy. This may be because CFD models are highly sensitive to local
geometric errors, while segmentation metrics are global and therefore may not
fully capture these localized deviations. Another reason may be that
differences in clipping (which were not accounted for by segmentation metrics)
are responsible for some of the CFD errors, as shown by our analysis of
equally clipped data. However, significant CFD errors remained after removing
this confounding factor, and these errors were still not correlated with
image-based segmentation metrics. Irrespective of the cause, the poor
correlation between segmentation and CFD errors has some important
implications. Specifically, in our application it might be better to combine
conventional global image-based losses with more CFD specific objective
measures during training.

Computational fluid dynamics can benefit in several ways from machine
learning. Firstly, ML segmentation is completely automated and very fast,
enabling significant reduction in pre-processing time, one of the major
impediments to clinical uptake. Secondly, ML segmentations are completely
reproducible, and this is important as we have shown significant human inter-
observer variability. Finally, there has been recent work demonstrating the
use of ML to accelerate the CFD simulations. Combined with ML segmentation
this would substantially reduce the time taken to perform CFD and make CFD
much more attractive for routine clinical use.

###  4.3 Limitations

One of the main limitations of this study was that a simplified CFD model was
applied across all subjects (laminar, steady state with no patient-specific
parameters). This was done to better isolate the effect of segmentation
differences on the resulting CFD model. However, it does limit the patient
specific aspect of these comparisons and in the future, boundary conditions
for each subject (such as velocity profiles taken from phase contrast MRI)
could be incorporated into the model. Furthermore, now that we have
demonstrated good agreement using simple CFD models, the utility of ML
segmentation for more complex CFD models should be investigated.

Another limitation is that the methods used for comparison don’t necessary
account for the full flow field. We used plane-averaged pressures and
velocities along the length of the centreline to quantitatively compare
different CFD models. However, this averaging does lead to a loss of localized
details in the flow fields. Additionally, the slice locations were determined
independently for ML, GT and SO models, so there may not be an exact one-to-
one correspondence. In future studies, particularly if using more complex CFD
models, new metrics of CFD errors that capture subtle deviations will need to
be developed.

##  5 Conclusions

A convolutional neural network was developed, optimized and trained for
segmentation of the aorta and the pulmonary arteries in 3D cardiovascular MRI.
The segmentation network was validated for its primary purpose: the creation
of CFD models and calculation of flow fields. Segmentation errors in terms of
Dice, IoU, HD and ASD as well as derived pressure and velocity field errors
were in the range of human inter-observer variability. The proposed method
could help to automate clinical hemodynamic assessment workflows and improve
their robustness.

##  6 Declarations

###  6.1 Ethics approval and consent to participate

The use of retrospectively collected training and test data was approved by
the local research ethics committee, and written consent was obtained from all
subjects/guardians (Ref: 06/Q0508/124).

###  6.2 Consent for publication

Written consent to use the data for publication was obtained from all
subjects/guardians (Ref: 06/Q0508/124).

###  6.3 Availability of data and materials

The datasets used and/or analysed during the current study are available from
the corresponding author on reasonable request.

###  6.4 Competing interests

The authors declare that they have no competing interests.

###  6.5 Funding

This work was supported by UK Research and Innovation (MR/S032290/1), Heart
Research UK (RG2661/17/20), the British Heart Foundation (NH/18/1/33511,
PG/17/6/32797), the Engineering and Physical Sciences Research Council
(EP/N02124X/1) and the European Research Council (ERC-2017-StG-757923).

This report incorporates independent research from the National Institute for
Health Research Biomedical Research Centre Funding Scheme. The views expressed
in this publication are those of the authors and not necessarily those of the
NHS, the National Institute for Health Research or the Department of Health.
The funders had no role in study design, collection, analysis and
interpretation of the data, decision to publish, or preparation of the
manuscript.

###  6.6 Authors’ contributions

  * •

JMT contributed to the conception and design of this work, pre-processed the
images and labels, designed, implemented, trained and optimized the neural
network and related functionality, extracted the reported metrics, interpreted
the results, performed the statistical analysis and was a major contributor in
drafting the manuscript.

  * •

EP contributed to the conception and design of this work, designed and
implemented the semi-automatic mesh processing pipeline, performed the
computational fluid dynamics simulations, extracted the reported pressure and
velocity field data, interpreted the results and was a major contributor in
drafting the manuscript.

  * •

RJ performed the reference ground truth segmentation.

  * •

ES contributed to the conception and design of this work.

  * •

RP contributed to the preparation of the external test set.

  * •

AS contributed to the preparation of the external test set.

  * •

CC contributed to the conception and design of this work.

  * •

JS contributed to the conception and design of this work and substantially
revised the manuscript.

  * •

SS contributed to the conception and design of this work and substantially
revised the manuscript.

  * •

VM contributed to the conception and design of this work, pre-processed the
images and labels, performed the secondary observer segmentation and
contributed to draft and substantially revise the manuscript.

All authors read and approved the final manuscript.

###  6.7 Acknowledgements

Not applicable.

##  Appendix A Sensitivity Analysis: Laminar vs Turbulent

It was found that an element count of 300,000 was well-suited to capture the
flow details in both the aorta and the PA, and for yielding a stable solution.
Iterations of 300 and 750 were sufficient for the aorta and PA respectively,
in order to reach convergence and accurate results. The error (MAPE) for the
aorta with 315,766 elements vs the highest resolution aorta with 1,102,454
elements was 2.02% in pressure, and 3.4% in velocity. The error (MAPE) for the
PA with 398,947 elements vs the highest resolution PA with 1,116,831 elements
was 3.23% in pressure, and 2.15% in velocity. Lastly, the Reynold’s number was
computed on both test cases, and resulted in 3221.8 and 4250.5 for the aorta
and pulmonary artery respectively. The same simulations were carried out with
a k-omega turbulence model (with default parameters), which showed no
difference when compared to the laminar model.

Figure 10: Difference velocity fields in Ao and PA when applying a k-omega
turbulence model(left = laminar, right = turbulent).

##  Appendix B Hyperparameter Optimization

Training, including hyperparameter optimization, took 24 h, during which a
total of 124 hyperparameter configurations were sampled by the Hyperband
iterations. The top 10 performing configurations in terms of mean validation
Dice score are reported in Table 2 and Fig. 11.

Table 2: Top 10 hyperparameter configurations. # | Scales |  Layers per block |  Initial filters | Learning rate |  Batch size | Loss function | Dice
---|---|---|---|---|---|---|---
1 | 3 |  2 |  64 | $3.46\cdot 10^{-4}$ |  2 | Focal Tversky | 0.946
2 | 3 |  4 |  32 | $1.99\cdot 10^{-4}$ |  2 | Tversky | 0.944
3 | 4 |  3 |  32 | $3.40\cdot 10^{-4}$ |  4 | Dice | 0.944
4 | 4 |  3 |  32 | $1.37\cdot 10^{-4}$ |  4 | IoU | 0.943
5 | 3 |  3 |  64 | $6.77\cdot 10^{-4}$ |  4 | IoU | 0.939
6 | 3 |  3 |  64 | $3.27\cdot 10^{-4}$ |  4 | Dice | 0.938
7 | 4 |  3 |  32 | $1.81\cdot 10^{-4}$ |  2 | Dice | 0.938
8 | 3 |  4 |  32 | $1.83\cdot 10^{-3}$ |  2 | Focal Tversky | 0.937
9 | 2 |  3 |  32 | $1.24\cdot 10^{-3}$ |  4 | Jaccard | 0.935
10 | 3 |  2 |  32 | $4.68\cdot 10^{-3}$ |  4 | Dice | 0.934
|  |  |  |  |  |  |
Figure 11: Parallel coordinates view of the top 10 hyperparameter
configurations. Each hyperparameter as well as the mean validation Dice score
is shown on its own axis. Each coloured line represents a hyperparameter
combination, with vertices at the corresponding values on the parallel axes.

The best performing configuration was as follows: scales = 3, layers per block
= 2, initial filters = 64, learning rate = $3.46\cdot 10^{-4}$, batch size =
2, and loss function = focal Tversky. This model was selected and used in all
further experiments.

The top 8 performing configurations have 3 or 4 scales, which suggests that
models with 2 scales may have too limited representational power or deep layer
receptive fields. Models with at least 3 layers per block also tended to
perform better, with the notable exception of the best performing model, which
had 2. A relatively small learning rate also seems to be advantageous, with 6
of the top 7 configurations below 0.0005. Confusion-based losses seem to
outperform cross-entropy, which does not appear on any of the top 10
configurations, but none of those appears to consistently outperform the rest.
Finally, although our best model had 64 initial filters, the abundance of top
performing models with 32 filters suggests it might be possible to find a
32-filter configuration with little or no performance penalty, which would
result in a 2x reduction in the number of parameters and computational cost.

##  Appendix C Manual and Equal Clipping

In all cases, the equally clipped data produces better agreement between the
ML and GT CFD simulations (Fig. 12). The median pressure and velocity error in
aortas is reduced through equal clipping by 2.1 and 1.0 percentage points,
respectively. The median pressure and velocity error in PAs is reduced through
equal clipping by 4.2 and 2.7 percentage points, respectively. The magnitude
of outliers is also reduced through equal clipping.

Figure 12: Pressure and velocity mean average percentage errors (MAPE) for
manual and equally clipped aortas and pulmonary arteries.

## References

  *   * (1) Biglino, G., Capelli, C., Bruse, J., Bosi, G.M., Taylor, A.M., Schievano, S.: Computational modelling for congenital heart disease: How far are we from clinical translation? BMJ Publishing Group (2017). https://doi.org/10.1136/heartjnl-2016-310423. http://heart.bmj.com/
  * (2) Wang, C., Frimmel, H., Smedby, O.: Fast level-set based image segmentation using coherent propagation. Medical physics 41(7) (2014). https://doi.org/10.1118/1.4881315
  * (3) Weese, J., Kaus, M., Lorenz, C., Lobregt, S., Truyen, R., Pekar, V.: Shape Constrained Deformable Models for 3D Medical Image Segmentation. In: Lecture Notes in Computer Science (including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) vol. 2082, pp. 380–387. Springer, ??? (2001). https://doi.org/10.1007/3-540-45729-1_38
  * (4) Peters, J., Ecabert, O., Meyer, C., Schramm, H., Kneser, R., Groth, A., Weese, J.: Automatic Whole Heart Segmentation in Static Magnetic Resonance Image Volumes. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 4792 LNCS(PART 2), 402–410 (2007). https://doi.org/10.1007/978-3-540-75759-7_49
  * (5) Chen, C., Qin, C., Qiu, H., Tarroni, G., Duan, J., Bai, W., Rueckert, D.: Deep Learning for Cardiac Image Segmentation: A Review. Frontiers in Cardiovascular Medicine 0, 25 (2020). https://doi.org/10.3389/FCVM.2020.00025
  * (6) Berhane, H., Scott, M., Elbaz, M., Jarvis, K., McCarthy, P., Carr, J., Malaisrie, C., Avery, R., Barker, A.J., Robinson, J.D., Rigsby, C.K., Markl, M.: Fully automated 3D aortic segmentation of 4D flow MRI for hemodynamic analysis using deep learning. Magnetic Resonance in Medicine 84(4), 2204–2218 (2020). https://doi.org/10.1002/MRM.28257
  * (7) Zhuang, X., Li, L., Payer, C., Štern, D., Urschler, M., Heinrich, M.P., Oster, J., Wang, C., Smedby, O., Bian, C., Yang, X., Heng, P.-A., Mortazi, A., Bagci, U., Yang, G., Sun, C., Galisot, G., Ramel, J.-Y., Brouard, T., Tong, Q., Si, W., Liao, X., Zeng, G., Shi, Z., Zheng, G., Wang, C., MacGillivray, T., Newby, D., Rhode, K., Ourselin, S., Mohiaddin, R., Keegan, J., Firmin, D., Yang, G.: Evaluation of algorithms for Multi-Modality Whole Heart Segmentation: An open-access grand challenge. Medical Image Analysis 58, 101537 (2019). https://doi.org/10.1016/J.MEDIA.2019.101537
  * (8) Bai, W., Sinclair, M., Tarroni, G., Oktay, O., Rajchl, M., Vaillant, G., Lee, A.M., Aung, N., Lukaschuk, E., Sanghvi, M.M., Zemrak, F., Fung, K., Paiva, J.M., Carapella, V., Kim, Y.J., Suzuki, H., Kainz, B., Matthews, P.M., Petersen, S.E., Piechnik, S.K., Neubauer, S., Glocker, B., Rueckert, D.: Automated cardiovascular magnetic resonance image analysis with fully convolutional networks. Journal of Cardiovascular Magnetic Resonance 2018 20:1 20(1), 1–12 (2018). https://doi.org/10.1186/S12968-018-0471-X
  * (9) Böttcher, B., Beller, E., Busse, A., Cantré, D., Yücel, S., Öner, A., Ince, H., Weber, M.-A., Meinel, F.G.: Fully automated quantification of left ventricular volumes and function in cardiac MRI: clinical evaluation of a deep learning-based algorithm. The International Journal of Cardiovascular Imaging 2020 36:11 36(11), 2239–2247 (2020). https://doi.org/10.1007/S10554-020-01935-0
  * (10) Steeden, J.A., Quail, M., Gotschy, A., Mortensen, K.H., Hauptmann, A., Arridge, S., Jones, R., Muthurangu, V.: Rapid whole-heart CMR with single volume super-resolution. Journal of Cardiovascular Magnetic Resonance 22(1), 56 (2020). https://doi.org/10.1186/s12968-020-00651-x
  * (11) Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional networks for biomedical image segmentation. In: Lecture Notes in Computer Science (including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 9351, pp. 234–241. Springer, ??? (2015). https://doi.org/10.1007/978-3-319-24574-4_28
  * (12) Li, L., Jamieson, K., Rostamizadeh, A., Talwalkar, A.: Hyperband: A novel bandit-based approach to hyperparameter optimization. Journal of Machine Learning Research 18(185), 1–52 (2018)
  * (13) Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., Kudlur, M., Levenberg, J., Monga, R., Moore, S., Murray, D.G., Steiner, B., Tucker, P., Vasudevan, V., Warden, P., Wicke, M., Yu, Y., Zheng, X.: Tensorflow: A system for large-scale machine learning. In: 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pp. 265–283. USENIX Association, Savannah, GA (2016). https://www.usenix.org/conference/osdi16/technical-sessions/presentation/abadi
  * (14) Montalt-Tordera, J.: TensorFlow MRI: A library of TensorFlow operators for computational magnetic resonance imaging. Zenodo (2022). https://doi.org/10.5281/zenodo.6897643
  * (15) He, K., Zhang, X., Ren, S., Sun, J.: Delving deep into rectifiers: Surpassing human-level performance on ImageNet classification. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 1026–1034. Institute of Electrical and Electronics Engineers Inc., ??? (2015). https://doi.org/10.1109/ICCV.2015.123
  * (16) Kingma, D.P., Ba, J.L.: Adam: A method for stochastic optimization. In: 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. International Conference on Learning Representations, ICLR, ??? (2015). https://arxiv.org/abs/1412.6980v9
  * (17) Silversmith, W.: cc3d: Connected Components on Multilabel 3D Images (2021). https://github.com/seung-lab/connected-components-3d/
  * (18) Gijsen, F.J.H., Van De Vosse, F.N., Janssen, J.D.: The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model. Journal of Biomechanics 32(6), 601–608 (1999). https://doi.org/10.1016/S0021-9290(99)00015-9
  * (19) Caballero, A.D., Laín, S.: A Review on Computational Fluid Dynamics Modelling in Human Thoracic Aorta. Cardiovascular Engineering and Technology 2013 4:2 4(2), 103–130 (2013). https://doi.org/10.1007/S13239-013-0146-6
  * (20) Bonfanti, M., Balabani, S., Greenwood, J.P., Puppala, S., Homer-Vanniasinkam, S., Díaz-Zuccarini, V.: Computational tools for clinical support: a multi-scale compliant model for haemodynamic simulations in an aortic dissection based on multi-modal imaging data. Journal of The Royal Society Interface 14(136) (2017). https://doi.org/10.1098/RSIF.2017.0632
  * (21) Gabe, I.T., Gault, J.H., Ross, J., Mason, D.T., Mills, C.J., Schillingford, J.P., Braunwald, E.: Measurement of instantaneous blood flow velocity and pressure in conscious man with a catheter-tip velocity probe. Circulation 40(5), 603–614 (1969). https://doi.org/10.1161/01.CIR.40.5.603
  * (22) Payer, C., Štern, D., Bischof, H., Urschler, M.: Multi-label Whole Heart Segmentation Using CNNs and Anatomical Label Configurations. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 10663 LNCS, 190–198 (2017). https://doi.org/10.1007/978-3-319-75541-0_20
  * (23) Sun, C., Shrivastava, A., Singh, S., Gupta, A.: Revisiting unreasonable effectiveness of data in deep learning era. In: Proceedings of the IEEE International Conference on Computer Vision, vol. 2017-October, pp. 843–852. Institute of Electrical and Electronics Engineers Inc., ??? (2017). https://doi.org/10.1109/ICCV.2017.97
  * (24) Krizhevsky, A., Sutskever, I., Hinton, G.E.: ImageNet Classification with Deep Convolutional Neural Networks. Advances in Neural Information Processing Systems 25 (2012)
  * (25) Bhalodia, R., Elhabian, S.Y., Kavan, L., Whitaker, R.T.: DeepSSM: A Deep Learning Framework for Statistical Shape Modeling from Raw Images. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 11167 LNCS, 244–257 (2018)
  * (26) Raju, A., Miao, S., Cheng, C.-T., Lu, L., Han, M., Xiao, J., Liao, C.-H., Huang, J., Harrison, A.P.: Deep Implicit Statistical Shape Models for 3D Medical Image Delineation (2021)
  * (27) Zhu, J.Y., Park, T., Isola, P., Efros, A.A.: Unpaired Image-to-Image Translation using Cycle-Consistent Adversarial Networks. Proceedings of the IEEE International Conference on Computer Vision 2017-October, 2242–2251 (2017). https://doi.org/10.48550/arxiv.1703.10593
  * (28) Yan, W., Wang, Y., Gu, S., Huang, L., Yan, F., Xia, L., Tao, Q.: The Domain Shift Problem of Medical Image Segmentation and Vendor-Adaptation by Unet-GAN (2019). https://doi.org/10.48550/arxiv.1910.13681
  * (29) Liu, L., Wolterink, J.M., Brune, C., Veldhuis, R.N.J.: Anatomy-aided deep learning for medical image segmentation: a review. Physics in Medicine & Biology 66(11), 11–01 (2021). https://doi.org/10.1088/1361-6560/ABFBF4