Two-dimensional materials and their heterostructures have opened up new possibilities for magnetism at the nanoscale. In this study, we utilize first-principles simulations to investigate the structural, electronic, and magnetic properties of Fe/WSe2/Pt systems containing pristine, defective, or doped WSe2 monolayers. The proximity effects of the ferromagnetic Fe layer are studied by considering defective and Vanadium-doped WSe2 monolayers. All heterostructures are found ferromagnetic, and the insertion of the transition metal dichalcogenide results in redistribution of spin orientation and an increased density of magnetic atoms due to the magnetized WSe2.  There is an increase in the overall total density of states at the Fermi level due to WSe2, however, the transition metal dichalcogenide may lose its distinct semiconducting properties due to the stronger than van der Waals coupling. Spin resolved electronic structure properties are linked to larger spin Seebeck coefficients found in heterostructures with WSe2 monolayers.

Heterostructures (HSTs) composed of ferromagnetic (FM) and heavy metal (HM) films are typical systems employing spin-charge conversion mechanisms for spintronics applications [Dieny et al., 2020]. The quality of such devices is controlled by the atomic composition of the interface and its spin-dependent properties. In this context, HSTs with 2D layered materials have become of great interest [Sierra et al., 2021]. Chemically inert layered materials can improve the interface atomic sharpness and they provide new ways of property tuning. For example, semiconducting transition metal dichalcogenides (TMDCs) have strong spin orbit coupling (SOC) which allows the manipulation of spins via electrical means [Xiao et al., 2012]. On the other hand, their strong spin photon coupling enables the control of magnetism with optical methods [Phan et al., 2021]. The recent discovery of FM ordering at room temperature in  WS2 and WSe2 monolayers doped by V atoms [Yun et al., 2020, Pham et al., 2020, Zhang et al., 2020, Ortiz Jimenez et al., 2021] has shown that TMDCs can be effective dilute 2D ferromagnets with robust magnetic states, which can also be light-tunable.

HSTs with semiconducting TMDC layers offer possibilities for new ways to change and control magnetic properties,  spin transport and spin conversion phenomena through novel interfacial coupling mechanisms. For example, much enhanced magnetic exchange fields, valley splitting, and polarization effects tunable via optical methods have been observed in WSe2/CrI3 HSTs [Zhong et al., 2017, Seyler et al., 2018, Zhong et al., 2020]. Several impressive experimental studies have shown that spin current generation via thermal gradients across FM/TMDC/HM structures is significantly altered when compared with FM/HM systems. Specifically, a large enhancement of the voltage difference due to the longitudinal spin Seebeck effect in Pt/Ni81Fe19 and Pt/YIG after inserting WS2 (WSe2) layers has been reported [Dastgeer et al., 2019, Lee et al., 2020, Lee et al., 2020, Lee et al., 2021]. It was also found that such an increase is dependent on the number of layers in the heterostructure, degree of coverage by the TMDC layers, and quality of the interface. Others have examined spin-dependent transport in Fe3O4/MoS2/Fe3O4 and Co/Al2O3/MoS2 junctions [Dankert et al., 2017, Galbiati et al., 2018]. Spin valves based on spin filtering in WS2/Co and NiFe/MoS2/NiFe have also been reported [Wang et al., 2015, Zatko et al., 2019]. A gate-tunable magnetic order in WSe2 by changing the occupation of the Se vacancy states are successfully demonstrated [Nguyen et al., 2021] and computational work investigates the origin of ferromagnetism in complex configurations of Se vacancies in V-doped WSe2 monolayers [Yun et al., 2022]. The generated Se vacancies play a key role in the structural change of PdSe2 because of the laser-induced phase transformation of layered PdSe2 to metallic phase PdSe_2–x [Shautsova et al., 2019]. Ozyilmaz et. al. suggested that the spin Hall effect in graphene/TMDC heterostructures may also originate from chalcogenide vacancies [Avsar et al., 2014].

In addition to these experimental demonstrations, the theoretical understanding of the interface properties and coupling mechanisms of FM/TMDC/HM systems must also be advanced. For example, the changes in the band structure  of individual FM, TMDC, and HM layers calculated from first principles can help one understand the overall behavior of the entire HST. Computed charge transfer, bonding, and interfacial hybridization are also necessary in order to highlight magnetic proximity effects unique to FM/TMDC/HM HSTs. Additionally, first principles simulations are necessary to reveal spin selection transfer mechanisms across doped and defective TMDCs.

In this paper, we present first principles simulations of Fe/WSe2/Pt HSTs to elucidate their basic electronic and magnetic properties for the first time. We consider several cases in which the TMDC is a pristine WSe2, a monolayer with V doping and a monolayer with various types of vacancies. The interplay between the proximity magnetism from the Fe substrate, the charge transfer and orbital hybdrization across the WSe2 monolayer and the large SOC effects from the Pt layer leads to a complex picture of interface magnetism. The calculated band structure, as well as the Fermi surface for each component of the considered HSTs, are instrumental for the detailed microscopic analysis of the electronic and spin-dependent behavior across the interface.

(a) Side view of Fe/Pt and Fe/WSe2/Pt HSTs. Top views of the various WSe2 monolayers as part of the Fe/WSe2/Pt HSTs: (b) pristine WSe2 (HST2) (c) WSe2(V_Se-Se), a defective monolayer with a pair of Se vacancies (HST3); (d) W_0.967V_0.033Se2, a defective monolayer with 3.3% V doping substituting W (HST4); (e) W_0.967V_0.033Se2(V_Se-Se), a defective monolayer with 3.3% V doping and a pair of Se vacancies (HST5); (f) WSe2(V_Se-Se + V_W), a defective monolayer with a pair of Se vacancies and a W vacancy (HST6). The interlayer distances are denoted in panel (a), while the electronic spin directions for all atoms of the monolayers are also shown with arrows upon completing the spin-polarized calculations.

Different HSTs have been constructed by taking the middle layer as pristine WSe2, as well as monolayers with V atom substitutional dopants or Se or W vacancies. A side view of the Fe/Pt HST (taken as a reference) and a HST with a WSe2 monolayer inserted in-between the Fe and Pt layers is shown in Fig. <ref>a.  As part of the Fe/WSe2/Pt HST, we consider a perfect WSe2 monolayer, a monolayer with two Se vacancies - WSe2 (V_Se-Se), a monolayer with 3.3% doping of V atoms substituting W atoms – W_0.957V_0.033Se2, a  monolayer with two Se vacancies and a 3.3% V doping - W_0.957V_0.033Se2 (V_Se-Se), and a monolayer containing two Se and one W vacancies – WSe2 (V_Se-Se + V_W), as shown in Fig. <ref>b-f.

In each case, the TMDC monolayer is in a 2H configuration. Fe and Pt layers with $(001)$ surface termination, have been simulated to create the HSTs having layer thickness 2.804 Å  and 3.926 Å, respectively. For the calculations, a supercell consisting of $6 \times 6$ Fe unit cells (108 atoms), $3\sqrt{3} \times 5$ unit cells for the perfect monolayer (90 atoms), and $3\sqrt{2} \times 3\sqrt{2}$ Pt unit cells (108 atoms) is constructed. In the case of the reference Fe/Pt HST, the TMDC monolayer is removed, and for the others, defective and/or doped WSe2 are part of each HST (as shown in Fig. <ref>).

Basic properties of the considered HSTs from Fig. <ref>. The strain $\epsilon$ (%) along the inequivalent $a, b$ lattice directions with and without the vdW correction is shown for the layer with the largest shrinkage.  The interlayer separations $d$, $d_1$, $d_2$ are also given with and without vdW interaction. The binding energy is calculated using the expression $E_b=\left(E_{HST}-\sum\limits_{i} E_i^{layer} \right)/N$, where $E_{HTS}$ is the total energy of the HST, $E_i^{layer}$ is the total energy of the $i^th$  isolated layer, and $N$ is the number of atoms in the HST unit cell. The Bader charge is shown for each layered component of the HSTs. The magnetization per atom $m_1$, $m_2$, $m_3$ corresponds to each layer and $M$ is the total average magnetization of each HST. The Bader charges and magnetizations are obtained with the vdW correction taken into the calculations.

2*     2c|Strain (%)                                                                                2c|Layer seperation (Å)                                                                                2c|$E_b$ (meV)             3c|Bader charge $(e)$                                                      1c|2*c]@c@$m_1$
($\mu_B$)     1c|2*c]@c@$m_2$
($\mu_B$)     1c|2*c]@c@$m_3$
($\mu_B$)     1c|2*c]@c@$M$

1c|No vdW                                                      1c|vdW                                    1c|No vdW                                                                  vdW                                                                     1c|No vdW       vdW          1c|Fe        1c|TMDC       1c|Pt     1c|                                                                               1c|                                                                               1c|
HST1                  1r|c]@r@ (a, Fe) -3.428
(b, Fe) -3.411     c]<EMAIL_ADDRESS>-3.620     1c|$d= 1.795$                                                              $d= 1.787$                                                              1c|-416.247     -509.738     1r|4.163     1c|–         -4.163                      2.579                                                                                1c|–                                                                             -0.058     2.521
HST2                  1r|c]@r@(a, Fe) -2.045
(b, Fe) -3.460     c]<EMAIL_ADDRESS>-3.800     1c|c]@c@$d_1= 2.405$
$d_2= 2.405$     c]@c@$d_1= 2.284$
$d_2= 2.289$     1c|-86.078      -195.852     1r|0.888     1r|0.888      -1.776                      2.588                                                                                0.014                                                                                0.023      2.625

HST3                  1r|c]@r@(a, Fe) -2.300
(b, Fe) -3.601     c]<EMAIL_ADDRESS>-3.906     1c|c]@c@$d_1= 2.363$
$d_2= 2.379$     c]@c@$d_1= 2.200$
$d_2= 2.279$     1c|-87.749      -197.653     1r|2.368     1r|-1.143     -1.226                      2.589                                                                                0.012                                                                                0.018     2.619
HST4                  1r|c]@r@(a. Fe) -1.994
(b, Fe) -3.392     c]<EMAIL_ADDRESS>-3.731     1c|c]@c@$d_1= 2.244$
$d_2= 2.320$     c]@c@$d_1= 2.143$
$d_2= 2.234$     1c|-88.237      -199.307     1r|3.450     1r|-2.796     -1.822                      2.589                                                                                0.026                                                                                0.022       2.637
HST5                  1r|c]@r@(a, Fe) -2.540
(b, Fe) -3.432     c]<EMAIL_ADDRESS>-3.777     1c|c]@c@$d_1= 2.311$
$d_2= 2.282$     c]@c@$d_1= 2.064$
$d_2= 2.209$     1c|-91.638      -201.194     1r|3.562     1r|-1.004     -2.559                      2.594                                                                                0.040                                                                                0.025      2.659
HST6                  1r|c]@r@(a, Fe) -2.255
(b, Fe) -3.358     c]<EMAIL_ADDRESS>-3.731     1c|c]@c@$d_1= 2.180$
$d_2= 2.191$     c]@c@$d_1= 2.043$
$d_2= 2.120$     1c|-98.094      -207.405     1r|4.196     1r|-1.001     -3.195                      2.596                                                                                0.016                                                                                0.025     2.637

The formation of the supercells introduces slight strains when compared with the lattice structures of the individual free-standing constituents. The strain $\epsilon$ is quantified using $\epsilon = \frac{a_{HST} - a_i}{a_i}  \times 100 \% \; $  ($a_{HST}$ - lattice constant of the HST; $a_i$ - lattice constant of the Fe, TMDC, or Pt layers). In Table <ref>, we show results for $\epsilon$ of the Fe layer along the inequivalent $a$ and $b$ lattice directions, and all data is given in Table S-1 in the Supporting Information. For the Fe and Pt layers, all values are negative indicating that the lattice has shrunk, and overall this  effect is more pronounced when the vdW interaction is taken into account. On the other hand, slight stretching along the $b$ direction and shrinkage along the $a$ direction are found for all considered $\ch{WSe2}$ monolayers. Nevertheless, $|\epsilon|$ does not exceed $ \sim 3.9\%$.

Interlayer separations for each HST  with and without the vdW corrections are also reported in Table <ref>. Including the vdW interaction in the simulations results in decreased $d_1$, $d_2$ as expected. In fact, the interlayer distances are smaller than the typical vdW separations (usually $\succeq 3$  Å), which is indicative of stronger interactions between the layers. Our results for the calculated binding energies (with included vdW correction), also given in Table <ref>, show that $E_b \sim -500$ meV for the reference Fe/Pt HST is increased to $E_b \sim -200$ meV for HSTs with a WSe2 monolayer, which is expected since the chemically inert  TMDC monolayer increases the separation between Fe and Pt. Nevertheless, $|E_b| \sim 200$ meV suggests that stronger chemical effects coexist with the vdW interaction typically associated with interlayer coupling of TMDCs [Pollard et al., 2016].

To elucidate the chemical interaction in the HSTs, the Bader charges (expressed as atomic charges in relation to the neutral ones) for each layer are also calculated and given in Table <ref>. The Bader charge is a reliable quantitative measure for the charge distribution in a given system [Henkelman et al., 2006, Koch et al., 2018]. From the numerical data in Table <ref>, we find that the Fe layer and the perfect WSe2 take equal amounts of positive charge from Pt. However, this balance is changed in the case of defective or doped WSe2 whose charge becomes negative, while the Fe layer acquires an enhanced positive charge. In all cases the Pt layer charge is always negative, but its magnitude is reduced when comparing with the Fe/Pt HST and the systems containing the TMDC layer. This is well correlated with the shorter $d$ and stronger $E_b$, as compared to $d_1$, $d_2$ and the binding energy for the HSTs involving a TMDC layer.

The electron gain or loss between the different components in the HSTs is intrinsically related to the atomic magnetic moment distribution, also given in Fig. <ref> and Table <ref>. The average magnetic moments $m_1, m_2, m_3$ after relaxation show how the spin polarization for each layer changes as different HSTs are constructed. In all cases, the calculated magnetic moments are along the direction perpendicular to the layer surfaces featuring FM polarization. For the reference Fe/Pt system, the overall magnetization is always positive, however, the Pt layer itself has opposite spin polarization compared to the one for Fe. In the Fe/WSe2/Pt (HST2), each atom from the TMDC monolayer is spin-polarized in the upward direction with an overall layer magnetization $0.014$ $\mu_B$ due to the proximity of the Fe layer. It is interesting to note that $m_1$ for the Fe layer has increased, while the Pt layer has reversed its spin polarization direction with a reduced by a factor of two magnitude when compared to the Fe/Pt HST. Nevertheless, the total magnetization in HSTs with WSe2 monolayers is enhanced when compared to the reference Fe/Pt. The effect is most pronounced for the case with the V doped defective monolayer due to all $m_1, m_2,$ and $m_3$ being the largest, as seen in  Table <ref>.

The spin direction of each atom (depicted in Fig. <ref>)   in the doped TMDC layer further shows that the V atom and its surrounding first and second nearest neighbors, even in the presence of Se vacancies, are coupled antiferromagnetically. This type of rather local AFM frustration is quenched, however, showing a preferred FM orientation over a longer period [Zhang et al., 2020]. On the other hand, all neighboring atoms in the presence of W and Se vacancies are FM coupled to their neighbors, as evident in HST6. The results indicate that the defective or doped HSTs have similar $m_1$, $m_2$, $m_3$, with the largest $m_2$ is observed for the W_0.967V_0.033Se2(V_Se-Se) from HST5 (Table <ref>).

An important factor for the bonding in the HSTs is the charge transfer, which can be calculated using $ \delta \rho_{HST} = \rho_{HST} - ( \rho_{Fe} + \rho_{WSe_2} + \rho_{Pt}) \text{ } $ whereas $\rho_{HST}$, $\rho_{Fe}$, $ \rho_{WSe_2}$, $\rho_{Pt}$ are the total charge in the HST, isolated Fe layer, isolated WSe2 monolayer and the isolated Pt layer, respectively. In Fig. <ref>, $\delta \rho_{HS}$  is shown for the considered HSTs from Fig. <ref>. Charge is moved from depleted regions (red isosurfaces) to regions with accumulated electrons (blue isosurfaces). A uniform layer of electron accumulation corresponds to ionic-like bonding as found at the Fe/Pt boundary indicating strong chemical interaction between the two layers (Fig. <ref>a). Upon insertion of the WSe2 monolayer, localized depletion above the Se atoms and localized electron accumulation above Pt and Fe atoms is observed (Fig. <ref>b). These ellipsoidal-like regions suggest dipole-dipole charge density distribution. The presence of Se vacancies (Fig. <ref>c) or V dopants (Fig. <ref>d) does not significantly affect the Fe/Se and Pt/Se interface as the main changes are primarily at the location of these disturbances. Comparing the different cases displayed in Fig. <ref> shows that the charge transfer is rather concentrated at the direct Fe/Pt interface, and the WSe2 monolayer serves as a medium through which charge is transferred weakens the bond between Fe and Pt layers.

Charge transfer in (a) Fe/Pt (HST1), (b) Fe/WSe2/Pt (HST2), (c) Fe/WSe2 (V_Se-Se)/Pt (HST3), (d) Fe/W_0.967V_0.033Se2/Pt (HST4), (e) Fe/W_0.967V_0.033Se2 (V_Se-Se)/Pt (HST5), (f) Fe/WSe2 (V_Se-Se + V_W)/Pt (HST6). The color coding for the atoms is the same as in Fig. <ref>. Total and atomically resolved density of states (DOS) for (g) HST1, (h) HST2, (i) HST3, (j) HST4, (k) HST5 and (l) HST6.

In Fig. <ref>g-l, the calculated density of states (DOS) is also shown. In the case of Fe/Pt, both types of atoms contribute to the DOS around the Fermi level. This situation does not change for the HSTs with WSe2 monolayers. However, we find that the presence of the TMDC always increases the total DOS at the Fermi level by about $23\%$ regardless of the presence of defects or dopants. This increase in conduction states results in a reduced conductivity mismatch between Fe and Pt, which can have a favorable effect on the enhancement of spin dependent transport phenomena in such systems, as observed experimentally [Dastgeer et al., 2019, Lee et al., 2020, Lee et al., 2020, Lee et al., 2021].

Further insight into the microscopic picture of such heterostructures can be gained by considering their energy band structures. Since the different layers are incommensurate, the constructed supercell contains several unit cells from each component. This situation also occurs in the case of monolayer WSe2 with impurities or V dopants, for which the supercell contains several unit cells. Because of the shrunk Brillouin zone (BZ) appearing in materials containing multiple unit cells, artificially folded bands are created. Thus, we are faced with the problem of band untangling on the primitive BZs for each layer to have a complete analysis of the electronic structure of the various HSTs. As representative examples, we consider Fe/Pt, Fe/WSe2/Pt and Fe/W_0.9V_0.1Se2/Pt to compare how the unpolarized perfect WSe2 and FM V-doped TMDC layer affect the band structure of the reference Fe/Pt.

Unfolded band structures for the: (a) Fe layer and (b) Pt layer from the Fe/Pt HST; (c) Fe layer, (d) Pt layer and (e) WSe2 monolayer from the Fe/WSe2/Pt HST; (f) standalone WSe2 monolayer; (g) Fe layer, (h) Pt layer, and (i) W_0.9V_0.1Se2 monolayer from the Fe/W_0.9V_0.1Se2/Pt HST; (j) standalone W_0.9V_0.1Se2 monolayer. The band structures are given in the native Brillouin zones for each material. The Brillouin zones (BZ) and high-symmetry points for the Fe, Pt and WSe2 layers are also shown.

In Fig. <ref>, we show the unfolded bands (obtained as discussed in Methods) as projected on each component native BZ along high symmetry points. For the Fe and Pt layers, the high symmetry points of $\Gamma$ (0 0 0) – X (0.5 0 0) – S (0.5 0.5 0) – Y (0 0 0.5) – $\Gamma$ are selected to calculate the energy bands. Fig. <ref>a, b show the orbitally projected bands for both layers. We find that the bands around the Fermi level along the $X - S - Y$ (Fig. <ref>a) and $\Gamma - X - S$ (Fig. <ref>b) paths are primarily composed of hybridized $d_{xy}$, $d_{yz}$, $d_{xz}$ orbitals. Along $\Gamma - X$ and $Y - \Gamma$ additional hybridization with $d_{z^2}$ and $d_{x^2-y^2}$ is also found. In the Supproting Information, a separate $p$; $d_{xy}$, $d_{yz}$, $d_{xz}$; $d_{z^2}$; $d_{x^2-y^2}$ band decomposition is given for easier viewing together with the unfolded band structures of the separate Fe/WSe_2 (Fig. S2) and Pt/WSe_2 (Fig. S3). Upon inserting WSe2, the orbital character of the bands around $E_F$ does not change, however, the Fe orbitals become less concentrated, while the contribution from the Pt orbitals and their intensity increase significantly (Fig. <ref>c, d, g, h).

It is further interesting to analyze the WSe2 monolayer: as a free-standing material, it is an indirect bandgap semiconductor with a gap of $1.390$ eV at the K – Q path, as shown in Fig. <ref>f, which is in good agreement with previous studies [Zhang et al., 2015, Hu et al., 2019]. The emergence of an indirect gap is a consequence of the SOC which is responsible for $d$-band splitting and lowering $d_{x^2-y^2}$ in the conduction region at the Q-point in Fig. <ref>f. The band structure of the monolayer changes significantly upon its insertion between Fe and Pt. Fig. <ref>e shows several bands crossing $E_F$ and a series of $d$-bands of various intensities around the K-point. It appears that the WSe2 monolayer has lost its semiconducting properties due to the stronger than vdW bonding with Fe and Pt, as also evidenced from Fig. <ref>e. The standalone doped monolayer W_0.9V_0.1Se2, on the other hand, shows an upward shift of the valence bands touching the Fermi level at the K-point. There is further splitting at the top of the valence maximum crossing over to the gap of the monolayer (Fig. <ref>j). Impurity levels originating from the V d-states are also found localized near the conduction band minimum. Nevertheless, the band structure of the W_0.9V_0.1Se2 as part of the HST in Fig. <ref>i washes away these signature features; instead, similar to WSe2 from Fe/WSe2/Pt  there is a series of bands near and at $E_F$ with various intensities  (Fig. <ref>e). The Fe and Pt layers, however, result in a similar situation as for the Fe/WSe2/Pt HST with a series of bands near and at $E_F$ with various intensities.

The unfolded Fermi Surface for the: (a) Fe layer and (b) Pt layer from the Fe/Pt heterostructure; (c) W_0.967V_0.033Se2 monolayer; (d) W_0.9V_0.1Se2 monolayer; (e) Fe layer, (f) Pt layer and (g) WSe2 monolayer from the Fe/WSe2/Pt heterostructure; (h) Fe layer, (i) Pt layer, and (j) W_0.9V_0.1Se2  from the Fe/W_0.9V_0.1Se2/Pt HST. The band structures are given in the native Brillouin zones for each material. The number of Fermi-surface states is given by the color bar.

To get a better understanding of the spin-resolved distribution across the interface, we also provide results for the computed Fermi Surfaces (FSs) for each HST component. The zone folding due to the reduced and incommensurate Brillouin zones in the supercell also appears in the FS for each HST. In order to monitor the evolution of the different carriers, the disentangled FSs for the majority and minority spins for each layer in Fe/Pt, Fe/WSe2/Pt and Fe/W_0.9V_0.1Se2/Pt HSTs in the $k_x-k_y$ plane are calculated. Fig. <ref>a shows rather separated regions in Fe designated to the majority and minority spins, while for Pt in Fig. <ref>b such a distinction is not much pronounced. Inserting the WSe2 layer changes the spin selection regions in $\mathbf{k}_{\perp}$ space for the Pt and Fe layers. Fig. <ref>f,i show that there are no minority-spin states close to the $X (0.5\; 0 \; 0)$, $Y (0 \; 0.5 \; 0)$ point whereas majority-spin states exist almost everywhere in the surface Brillouin zone.

We next examine the FS of the TMDC monolayer. Fig. <ref>c, d shows that extended regions (around $K (1/3 \; 1/3 \; 0)$ points) are designated for both majority carriers despite the fact, that an isolated WSe2 does not have an inherent magnetism. The isolated V-doped monolayer, on the other hand, displays regions centered around the $K$-point for majority spins only, whose spread increases as the V-doping is increased (Fig. <ref>c, d). The FS for the W_0.9V_0.1Se2 from the Fe/W_0.9V_0.1Se2/Pt becomes very similar to the one for WSe2 from Fe/WSe2/Pt, which shows that the spin transfer mechanism through the TMDC does not change significantly upon doping.

The spin dependent transport properties of the HSTs are directly related to their electronic structure. The spin Seebeck effect associated with the thermal injection of spin currents is of particular interest for spintronics and spin caloritronic devices, which rely on large spin Seebeck coefficients[Adachi et al., 2013].  Using linear response theory within the transport Boltzmann equation, we calculate the spin Seebeck coefficient
($\sigma_{maj,min}, S_{maj,min}$ - electric conductivities and Seebeck coefficients for the majority and minority carriers) for several HSTs. Fig. <ref>a shows that the magnitude of $S_{spin}$ at $T=300$ K for the HSTs with WSe2 is larger than the one for Fe/Pt for practically all values of the chemical potential. For example, at $\mu=0$ $|S_{spin}|$ for Fe/W_0.9V_0.1Se2/Pt is more than doubled the one for Fe/Pt. In the peak-like features in the spin Seebeck coefficients at $\mu\sim -0.1$ eV and $\mu\sim -0.8$ eV, this increase is even strong as $|S_{spin}|$ for Fe/W_0.9V_0.1Se2/Pt shows triple values compared to Fe/Pt.

One of the factors contributing to the overall enhancement of the spin Seebeck coefficient magnitude is the decreased $\sigma_{maj}+\sigma_{min}$ for the HSTs with the TMDC monolayer as seen in Fig. <ref>e-g. The insertion of WSe2 leads to reduced charge transfer due to the larger separation between the Fe and Pt layers (Table <ref>), which explains the decrease in the overall electronic transport for Fe/WSe2/Pt and Fe/W_0.9V_0.1Se2/Pt. Fig. <ref>b-d further shows that there is larger disparity between $S_{maj}$ and $S_{min}$ for the HSTs with the TMDC as compared with Fe/Pt for which the spin resolved Seebeck coefficients have similar magnitudes for most $\mu$. Our calculations also show that the enhanced difference between $|S_{maj}|$ and $|S_{min}|$ is attributed to the larger spin majority channels for HSTs with WSe2 as found from the FS distributions in Fig. <ref>. Given that Seebeck coefficients are associated with localized features in DOS [Popescu and Woods, 2010, Mahan and Sofo, 1996], the larger peaks in the majority DOS (Fig. <ref>) is another factor contributing to the enhanced difference in the spin resolved Seebeck coefficients for the HSTs with TMDC monolayers.

(a) The spin Seebeck coefficient $S_{spin}$ as a function of chemical potential $\mu$ for the Fe/Pt, Fe/WSe_2/Pt, Fe/W_0.9V_0.1Se2/Pt HSTs;   The majority $S_{maj}$ and minority $S_{min}$ Seebeck coefficients as a function of $\mu$ for (b) Fe/Pt, (c) Fe/WSe2/Pt, (d) Fe/W_0.9V_0.1Se2/Pt HSTs; The majority $\sigma_{maj}$, minority
$\sigma_{min}$ and total $\sigma_{maj}+\sigma_{min}$ electric conductivities divided by the constant relaxation time $\tau_0$ as a function of
$\mu$ for (b) Fe/Pt, (c) Fe/WSe2/Pt, (d) Fe/W_0.9V_0.1Se2/Pt HSTs. All results are at $T=300$ K.

Our first principles simulations show that  Fe/WSe2/Pt HSTs exhibit complex electronic and magnetic properties. The pristine WSe2, which is non magnetic, acquires magnetic moments on its atoms due to its proximity to the FM substrate of Fe. This distribution can be locally disturbed when defects or dopants are present. The chemically inert WSe2 weakens the bond between Fe and Pt, but significant charge transfer occurs through the monolayer in the HST. This rather stronger than vdW coupling may be seen in the DOS and band structure results showing the detailed orbital hybridization balance. The insertion of the TMDC monolayer also leads to an increased DOS at the Fermi level, which is about   $(23\%)$ higher than DOS of the reference HST. This is interesting to compare with  YIG/C60/Pt: a related HST in which the C60 layer is also a semiconductor which bonds relatively strongly to YIG and Pt [Kalappattil et al., 2020]. In this case, the six-fold DOS increase at the Fermi level is attributed to a much reduced conductivity mismatch, which explained the giant longitudinal spin Seebeck in YIG/C60/Pt when compared with YIG/Pt [Kalappattil et al., 2020]. DOS increase in Fe/Pt with TMDC monolayers is also present, which can be a factor responsible for the longitudinal spin Seebeck enhancement as reported in [Dastgeer et al., 2019, Lee et al., 2020, Lee et al., 2020, Lee et al., 2021].

The unfolded band structures for the Fe/WSe2/Pt HSTs are especially useful in understanding to what extent each component preserves its electronic properties. Of special interest is the WSe2 and what happens to its semiconducting nature after its insertion between the Fe and Pt substrates. Our calculations show that much of the “identity” of the TMDC monolayer is lost as there are numerous bands with various intensities crossing the Fermi level. The strong hybridization with orbitals from the substrates alters the character of the energy bands of the  WSe2, which shows that its semiconducting behavior has changed. The presence of defects or dopants do not alter significantly this type of overall behavior. In light of this situation, potential experiments accessing the energy gap of the TMDC from the HST via optical means may need to be reevaluated [Phan et al., 2021].

The spin resolved atomic distribution and the FS calculations are also instrumental in understanding how localized spin is transferred across the interface of Fe/WSe2/Pt HSTs. Due to the proximity of the Fe substrate,  WSe2 becomes ferromagnetically polarized. The Pt substrate also has magnetization $m_3\neq 0$. Furthermore, the overall magnetization of Fe/WSe2/Pt is enhanced when comparing with the reference Fe/Pt HST. The spin resolved FSs show large overlapping regions in $\mathbf{k_i}$ space for both spin carriers, and for the WSe2 there is little difference between the areas corresponding to the majority or minority carriers. As a result, we see that there is an increased density of magnetic atoms giving rise to a higher $M$, which effectively enhances the spin mixing and overall spin transfer in the HST as shown theoretically in [Jia et al., 2011].

Our first principles simulations of Fe/WSe2/Pt with pristine, defective or V-doped TMDC monolayer reveal complex electronic and magnetic properties of the HSTs and their individual components. Although the ferromagnetism in WSe2 can be tuned via doping and/or defects, we find that this effect becomes secondary in the overall behavior of the HST. The semiconducting nature of the TMDC becomes diluted when the monolayer is between the Fe and Pt substrates. The total DOS increase at the Fermi level, larger magnetization, and importantly much enhanced spin density at the interface with overlapping majority and minority regions in the FS of Fe/WSe2/Pt HST influence the electronic and spin transport across the interface. The combined effect of all of these factors results in larger spin Seebeck coefficients in HSTs with TMDC layers, which can be useful in interpreting recent experiments reporting enhanced magnetization and spin resolved transport properties [Yun et al., 2020, Pham et al., 2020, Zhang et al., 2020, Dastgeer et al., 2019, Lee et al., 2020, Lee et al., 2020, Lee et al., 2021, Wang et al., 2015, Nguyen et al., 2021, Yun et al., 2022].

Computational Methods: The first principles simulations are performed with the DFT Vienna ab initio Simulations Package (VASP), which is based on projector-augmented wave (PAW) pseudopotentials and periodic boundary conditions [Kresse and Furthmüller, 1996, Kresse and Furthmüller, 1996]. The exchange-correlation effects are considered using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation with an included spin-orbit coupling (SOC) [Perdew et al., 1996]. The vdW interaction is included via the vdW-D3 (BJ) dispersion correction [Grimme et al., 2011]. The electronic shell structure is taken to have 8 valence electrons of Fe (3d^6 4s^2), 6 valence electrons of W (6s^2 5d^4), 6 valence electrons of Se (4s^2 4p^4) and 10 valence electrons of Pt (5d^9 6s^1). The kinetic energy cutoff is $450$ eV. In the construction of the HST supercells, we use the rectangular unit cell of WSe_2 $(a = 3.32 \text{ \AA}, b = 5.75 \text{ \AA})$ with a $6 \times 6$ Fe supercell (108 Fe atoms), a $3 \sqrt{3} \times 5$ WSe_2 supercell ($30$ W and $60$ Se atoms), a $3 \sqrt{2} \times 3 \sqrt{2}$ Pt supercell ($108$ Pt atoms). To eliminate the interactions between periodic images, a vacuum layer of $15$ Å is taken between the HST layers for all calculations.
The relaxation criteria of all systems is performed by allowing all cell parameters to change with $10^{-5}$  eV energy and $0.01$ eV⁄Å force relaxation criteria. The intergration over the Brillouin zone is performed at a single $\Gamma$ point for the optimization calculation. The DOS, charge density differences and Bader charges calculations are accomplished with $2 \times 2 \times 1$ $\Gamma$-centered Monkhorst-Pack k-point mesh. These calculations are performed with employing the Gaussian smearing with 0.1 eV width. To account for the strongly correlated electronic systems, we employed the $GGA+U$ method within Dudarev’s approach with an on-site Coulomb correction in the localized d orbitals of the V atoms. The parameter value of $U_{eff} = U - J = 3.0$ eV is selected, which has been shown to describe the ground-state electronic and magnetic properties of the system accurately [Duong et al., 2019, Yun et al., 2020]. For the Fe, W, Se and Pt, $U_{eff} = 0$ eV. The visualization of all systems is attained in VESTA [Momma and Izumi, 2011].

Given the zone folding problem that occurs in simulations involving supercells, the disentanglement of the various band structures is performed by adapting the vaspkit code [Wang et al., 2021]. Within this code, one creates $\mathbf{k_i}$ vectors along the high symmetry directions in the primitive Brillouin zone and then translates them into the supercell reciprocal space. Performing VASP calculation using these $\mathbf{k_i}$ vectors yields the plane wave coefficients and eigenvalues of each state necessary to determine spectral weights for the projection of the supercell eigenstates on the primitive cell eigenstates.

The FS unfolding is achieved using the bands4vasp package [Dirnberger et al., 2021], in which  the energy eigenvalues for every $\mathbf{k_i}$ primitive cell vector are assigned to a particular energy band by considering similarities in the orbital symmetry of the eigenstates. Then, a band center of mass is defined by averaging the energy values with weights given by the corresponding $P_{ \mathbf{K} m} ( \mathbf{k_i})$ character corresponding to the  projection of the supercell eigenstates $\left| \mathbf{K} m \right>$ on the primitive cell eigenstate $\left| \mathbf{k} n \right>$. Finally, the intersection of every band with the Fermi level is found by interpolation. Because the band structure and FS unfolding are computationally expensive, we consider the supercells consisting of $2 \times 6$ Fe unit cells (36 atoms), $\sqrt{3} \times 5$ unit cells for the  monolayer (30 atoms), and $\sqrt{2} \times 3\sqrt{2}$ Pt unit cells (36 atoms) are constructed. The optimization calculations of these system are performed with $6 \times 2 \times 1$ $\Gamma$-centered Monkhorst-Pack k-point mesh.

The transport properties of the HSTs are calculated with the BoltzTraP2 code [Madsen et al., 2018] with a large k-mesh of 112 000 points for more accurate results. The code utilizes linear response theory of the Boltzmann equation within a constant relaxation time approximation. In general, this is a good approximation as long as the energy changes are at a smaller scale as compared to the energy changes in the density of states, which is typically the case for metals and many semiconductors.

D. T.-X. D. would like to acknowledge support from the Presidential Fellowship sponsored by the University of South Florida. L.M.W. acknowledges financial support from the US Department of Energy under Grant No.DE-FG02-06ER46297. Computational resources are provided by USF Research Computing.

The Supporting Information contains structural data for the strains in the considered heterostructures upon the constructuion their supercells. Results for the calculated unfolded band structures of Fe/Pt, Fe/WSe2, Pt/WSe2, Fe/WSe2/Pt and Fe/W_0.9V_0.1Se2/Pt heterostructures and their constituents are also shown for further comparison with the results presented in the main text.

The Supporting Information is available free of charge at <https://pubs.acs.org/doi/10.1021/acs.jpclett.2c01925>.

[Dieny et al., 2020]
Dieny, B.; Prejbeanu, I. L.; Garello, K.; Gambardella, P.; Freitas, P.;
Lehndorff, R.; Raberg, W.; Ebels, U.; Demokritov, S. O.; Akerman, J.,
et al.  Opportunities and challenges for spintronics in the
microelectronics industry. Nature Electronics 2020, 3,
[Sierra et al., 2021]
Sierra, J. F.; Fabian, J.; Kawakami, R. K.; Roche, S.; Valenzuela, S. O. Van
der Waals heterostructures for spintronics and opto-spintronics. Nature
Nanotechnology 2021, 16, 856–868
[Xiao et al., 2012]
Xiao, D.; Liu, G.-B.; Feng, W.; Xu, X.; Yao, W. Coupled Spin and Valley Physics
in Monolayers of ${\mathrm{MoS}}_{2}$ and Other Group-VI Dichalcogenides.
Physical Review Letters 2012, 108, 196802
[Phan et al., 2021]
Phan, M.-H.; Trinh, M. T.; Eggers, T.; Kalappattil, V.; Uchida, K.-i.;
Woods, L. M.; Terrones, M. A perspective on two-dimensional van der Waals
opto-spin-caloritronics. Applied Physics Letters 2021,
119, 250501
[Yun et al., 2020]
Yun, S. J.; Duong, D. L.; Ha, D. M.; Singh, K.; Phan, T. L.; Choi, W.;
Kim, Y.-M.; Lee, Y. H. Ferromagnetic Order at Room Temperature in Monolayer
WSe2 Semiconductor via Vanadium Dopant. Advanced Science
2020, 7, 1903076
[Pham et al., 2020]
Pham, Y. T. H.; Liu, M.; Jimenez, V. O.; Yu, Z.; Kalappattil, V.; Zhang, F.;
Wang, K.; Williams, T.; Terrones, M.; Phan, M.-H. Tunable Ferromagnetism and
Thermally Induced Spin Flip in Vanadium-Doped Tungsten Diselenide Monolayers
at Room Temperature. Advanced Materials 2020, 32,
[Zhang et al., 2020]
Zhang, F. et al.  Monolayer Vanadium-Doped Tungsten Disulfide: A
Room-Temperature Dilute Magnetic Semiconductor. Advanced Science
2020, 7
[Ortiz Jimenez et al., 2021]
Ortiz Jimenez, V.; Pham, Y. T. H.; Liu, M.; Zhang, F.; Yu, Z.; Kalappattil, V.;
Muchharla, B.; Eggers, T.; Duong, D. L.; Terrones, M.; Phan, M.-H.
Light-Controlled Room Temperature Ferromagnetism in Vanadium-Doped Tungsten
Disulfide Semiconducting Monolayers. Advanced Electronic Materials
2021, 7, 2100030
[Zhong et al., 2017]
Zhong, D.; Seyler, K. L.; Linpeng, X.; Cheng, R.; Sivadas, N.; Huang, B.;
Schmidgall, E.; Taniguchi, T.; Watanabe, K.; McGuire, M. A., et al.
Van der Waals engineering of ferromagnetic semiconductor heterostructures for
spin and valleytronics. Science advances 2017, 3,
[Seyler et al., 2018]
Seyler, K. L.; Zhong, D.; Huang, B.; Linpeng, X.; Wilson, N. P.; Taniguchi, T.;
Watanabe, K.; Yao, W.; Xiao, D.; McGuire, M. A.; Fu, K.-M. C.; Xu, X. Valley
Manipulation by Optically Tuning the Magnetic Proximity Effect in
WSe2/CrI3 Heterostructures. Nano Letters 2018,
18, 3823–3828
[Zhong et al., 2020]
Zhong, D.; Seyler, K. L.; Linpeng, X.; Wilson, N. P.; Taniguchi, T.;
Watanabe, K.; McGuire, M. A.; Fu, K.-M. C.; Xiao, D.; Yao, W., et al.
Layer-resolved magnetic proximity effect in van der Waals heterostructures.
Nature Nanotechnology 2020, 15, 187–191
[Dastgeer et al., 2019]
Dastgeer, G.; Shehzad, M. A.; Eom, J. Distinct Detection of Thermally Induced
Spin Voltage in Pt/WS2/Ni81Fe19 by the Inverse Spin Hall Effect.
ACS Applied Materials & Interfaces 2019, 11,
[Lee et al., 2020]
Lee, S.-K.; Lee, W.-Y.; Kikkawa, T.; Le, C. T.; Kang, M.-S.; Kim, G.-S.;
Nguyen, A. D.; Kim, Y. S.; Park, N.-W.; Saitoh, E. Enhanced Spin Seebeck
Effect in Monolayer Tungsten Diselenide Due to Strong Spin Current Injection
at Interface. Advanced Functional Materials 2020, 30,
[Lee et al., 2020]
Lee, W.-Y.; Park, N.-W.; Kang, M.-S.; Kim, G.-S.; Jang, H. W.; Saitoh, E.;
Lee, S.-K. Surface Coverage Dependence of Spin-to-Charge Current across
Pt/MoS2/Y3Fe5O12 Layers via Longitudinal Spin Seebeck Effect. The
Journal of Physical Chemistry Letters 2020, 11,
[Lee et al., 2021]
Lee, W.-Y.; Kang, M.-S.; Kim, G.-S.; Park, N.-W.; Choi, K.-Y.; Le, C. T.;
Rashid, M. U.; Saitoh, E.; Kim, Y. S.; Lee, S.-K. Role of Ferromagnetic
Monolayer WSe2 Flakes in the Pt/Y3Fe5O12 Bilayer Structure in the
Longitudinal Spin Seebeck Effect. ACS Applied Materials & Interfaces
2021, 13, 15783–15790
[Dankert et al., 2017]
Dankert, A.; Pashaei, P.; Kamalakar, M. V.; Gaur, A. P. S.; Sahoo, S.;
Rungger, I.; Narayan, A.; Dolui, K.; Hoque, M. A.; Patel, R. S.;
de Jong, M. P.; Katiyar, R. S.; Sanvito, S.; Dash, S. P. Spin-Polarized
Tunneling through Chemical Vapor Deposited Multilayer Molybdenum Disulfide.
ACS Nano 2017, 11, 6389–6395
[Galbiati et al., 2018]
Galbiati, M.; Vecchiola, A.; Mañas-Valero, S.; Canet-Ferrer, J.; Galceran, R.;
Piquemal-Banci, M.; Godel, F.; Forment-Aliaga, A.; Dlubak, B.; Seneor, P.;
Coronado, E. A Local Study of the Transport Mechanisms in MoS2 Layers
for Magnetic Tunnel Junctions. ACS Applied Materials & Interfaces
2018, 10, 30017–30021
[Wang et al., 2015]
Wang, W.; Narayan, A.; Tang, L.; Dolui, K.; Liu, Y.; Yuan, X.; Jin, Y.; Wu, Y.;
Rungger, I.; Sanvito, S.; Xiu, F. Spin-Valve Effect in NiFe/MoS2/NiFe
Junctions. Nano Letters 2015, 15, 5261–5267
[Zatko et al., 2019]
Zatko, V. et al.  Band-Structure Spin-Filtering in Vertical Spin Valves
Based on Chemical Vapor Deposited WS2. ACS Nano 2019,
13, 14468–14476
[Nguyen et al., 2021]
Nguyen, T. D.; Jiang, J.; Song, B.; Tran, M. D.; Choi, W.; Kim, J. H.;
Kim, Y.-M.; Duong, D. L.; Lee, Y. H. Gate-Tunable Magnetism via Resonant
Se-Vacancy Levels in WSe2. Advanced Science 2021,
8, 2102911
[Yun et al., 2022]
Yun, S. J.; Cho, B. W.; Dinesh, T.; Yang, D. H.; Kim, Y. I.; Jin, J. W.;
Yang, S.-H.; Nguyen, T. D.; Kim, Y.-M.; Kim, K. K.; Duong, D. L.; Kim, S.-G.;
Lee, Y. H. Escalating Ferromagnetic Order via Se-Vacancies Near Vanadium in
WSe2 Monolayers. Advanced Materials 2022, 34,
[Shautsova et al., 2019]
Shautsova, V.; Sinha, S.; Hou, L.; Zhang, Q.; Tweedie, M.; Lu, Y.; Sheng, Y.;
Porter, B. F.; Bhaskaran, H.; Warner, J. H. Direct Laser Patterning and Phase
Transformation of 2D PdSe2 Films for On-Demand Device Fabrication.
ACS Nano 2019, 13, 14162–14171
[Avsar et al., 2014]
Avsar, A.; Tan, J. Y.; Taychatanapat, T.; Balakrishnan, J.; Koon, G.; Yeo, Y.;
Lahiri, J.; Carvalho, A.; Rodin, A.; O’Farrell, E., et al.
Spin–orbit proximity effect in graphene. Nature communications
2014, 5, 1–6
[Pollard et al., 2016]
Pollard, T. D.; Earnshaw, W. C.; Lippincott-Schwartz, J.; Johnson, G.
Cell biology E-book; Elsevier Health Sciences, 2016
[Henkelman et al., 2006]
Henkelman, G.; Arnaldsson, A.; Jónsson, H. A fast and robust algorithm for
Bader decomposition of charge density. Computational Materials Science
2006, 36, 354–360
[Koch et al., 2018]
Koch, D.; Golub, P.; Manzhos, S. Stability of charges in titanium compounds and
charge transfer to oxygen in titanium dioxide. Journal of Physics:
Conference Series 2018, 1136, 012017
[Zhang et al., 2015]
Zhang, C.; Chen, Y.; Johnson, A.; Li, M.-Y.; Li, L.-J.; Mende, P. C.;
Feenstra, R. M.; Shih, C.-K. Probing Critical Point Energies of Transition
Metal Dichalcogenides: Surprising Indirect Gap of Single Layer WSe2.
Nano Letters 2015, 15, 6494–6500
[Hu et al., 2019]
Hu, F.; Tao, L.; Ye, H.; Li, X.; Chen, X. ZnO/WSe2 vdW heterostructure for
photocatalytic water splitting. Journal Materials Chemistry C
2019, 7, 7104–7113
[Adachi et al., 2013]
Adachi, H.; Uchida, K.-i.; Saitoh, E.; Maekawa, S. Theory of the spin Seebeck
effect. Reports on Progress in Physics 2013, 76,
[Popescu and Woods, 2010]
Popescu, A.; Woods, L. M. Enhanced thermoelectricity in composites by
electronic structure modifications and nanostructuring. Applied Physics
Letters 2010, 97, 052102
[Mahan and Sofo, 1996]
Mahan, G. D.; Sofo, J. O. The best thermoelectric. Proceedings of the
National Academy of Sciences 1996, 93, 7436–7439
[Kalappattil et al., 2020]
Kalappattil, V.; Geng, R.; Das, R.; Pham, M.; Luong, H.; Nguyen, T.;
Popescu, A.; Woods, L. M.; Kläui, M.; Srikanth, H.; Phan, M. H. Giant spin
Seebeck effect through an interface organic semiconductor. Materials
Horizons 2020, 7, 1413–1420
[Jia et al., 2011]
Jia, X.; Liu, K.; Xia, K.; Bauer, G. E. W. Spin transfer torque on magnetic
insulators. Europhysics Letters 2011, 96, 17005
[Kresse and Furthmüller, 1996]
Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio
total-energy calculations using a plane-wave basis set. Physical Review
B 1996, 54, 11169–11186
[Kresse and Furthmüller, 1996]
Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations
for metals and semiconductors using a plane-wave basis set.
Computational Materials Science 1996, 6, 15–50
[Perdew et al., 1996]
Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made
Simple. Physical Review Letter 1996, 77,
[Grimme et al., 2011]
Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in
dispersion corrected density functional theory. Journal of
Computational Chemistry 2011, 32, 1456–1465
[Duong et al., 2019]
Duong, D. L.; Yun, S. J.; Kim, Y.; Kim, S.-G.; Lee, Y. H. Long-range
ferromagnetic ordering in vanadium-doped WSe2 semiconductor.
Applied Physics Letters 2019, 115, 242406
[Momma and Izumi, 2011]
Momma, K.; Izumi, F. VESTA3 for three-dimensional visualization of
crystal, volumetric and morphology data. Journal of Applied
Crystallography 2011, 44, 1272–1276
[Wang et al., 2021]
Wang, V.; Xu, N.; Liu, J.-C.; Tang, G.; Geng, W.-T. VASPKIT: A user-friendly
interface facilitating high-throughput computing and analysis using VASP
code. Computer Physics Communications 2021, 267,
[Dirnberger et al., 2021]
Dirnberger, D.; Kresse, G.; Franchini, C.; Reticcioli, M. Electronic State
Unfolding for Plane Waves: Energy Bands, Fermi Surfaces, and Spectral
Functions. The Journal of Physical Chemistry C 2021,
125, 12921–12928
[Madsen et al., 2018]
Madsen, G. K.; Carrete, J.; Verstraete, M. J. BoltzTraP2, a program for
interpolating band structures and calculating semi-classical transport
coefficients. Computer Physics Communications 2018,
231, 140–145


# Robotic Sonographer: Autonomous Robotic Ultrasound using Domain Expertise in
Bayesian Optimization

Deepak Raina12∗, SH Chandrashekhara3, Richard Voyles2, Juan Wachs2, Subir
Kumar Saha1
This work was supported in part by SERB (India) - OVDF Award No.
SB/S9/Z-03/2017-VIII; PMRF - IIT Delhi under Ref. F.No.35-5/2017-TS.I:PMRF;
National Science Foundation (NSF) USA under Grant #2140612; Daniel C. Lewis
Professorship and PU-IUPUI Seed Grant.1Indian Institute of Technology (IIT),
Delhi, India ({deepak.raina, saha}@mech.iitd.ac.in); 2Purdue University (PU),
Indiana, USA ({draina, rvoyles, jpwachs}@purdue.edu); 3All India Institute of
Medical Sciences (AIIMS), Delhi, India
(drchandruradioaiims@gmail.com).∗Corresponding author is Deepak Raina

###### Abstract

Ultrasound is a vital imaging modality utilized for a variety of diagnostic
and interventional procedures. However, an expert sonographer is required to
make accurate maneuvers of the probe over the human body while making sense of
the ultrasound images for diagnostic purposes. This procedure requires a
substantial amount of training and up to a few years of experience. In this
paper, we propose an autonomous robotic ultrasound system that uses Bayesian
Optimization (BO) in combination with the domain expertise to predict and
effectively scan the regions where diagnostic quality ultrasound images can be
acquired. The quality map, which is a distribution of image quality in a
scanning region, is estimated using Gaussian process in BO. This relies on a
prior quality map modeled using expert’s demonstration of the high-quality
probing maneuvers. The ultrasound image quality feedback is provided to BO,
which is estimated using a deep convolution neural network model. This model
was previously trained on database of images labelled for diagnostic quality
by expert radiologists. Experiments on three different urinary bladder
phantoms validated that the proposed autonomous ultrasound system can acquire
ultrasound images for diagnostic purposes with a probing position and force
accuracy of $\mathbf{98.7\%}$ and $\mathbf{97.8\%}$, respectively.

##  I Introduction

Ultrasound is the most frequently used imaging modality for diagnostic and
surgical interventions due to its low cost, non-ionizing nature, portability
and real-time feedback. Ultrasound offers several advantages over other
imaging modalities, like Magnetic Resonance Imaging (MRI) and Computed
Tomography (CT), however, the diagnosis by ultrasound is a highly operator-
dependent modality [1]. This is because of the skills required for manual
control of the probe and quality assessment of acquired images. Sonographers
employ both directed as well as random explorations strategies to search for
diagnostic-quality images. The ultrasound probe is moved within the region of
interest through hand maneuvers initially and fine adjustments to the probe’s
translational and rotational motion later. These maneuvers also include the
safe and precise adjustment of the pressure through the probe while
simultaneously analyzing the quality of acquired images. Such an intricate
procedure requires a great deal of skill, focus, experience and manual effort
from sonographers. In rural settings, skilled sonographers availability is
limited [2], and alternative solutions are required.

In order to reduce the burden on experts, a Robotic Ultrasound System (RUS) is
introduced. RUS consists of a dexterous robotic arm and an ultrasound machine
with its probe attached to the end effector of the robot, as shown in Fig. 1.
RUS can help ensure the accuracy, safety and consistency of the ultrasound
procedures. Recently, in order to address the aforementioned needs, several
telerobotic or human-assisted ultrasound systems have been proposed [3, 4, 5,
6, 7].

Figure 1: Robotic ultrasound system with probe attached to its end-effector
[3], conducting a urinary bladder ultrasound.

Compared to these systems, a fully automated ultrasound system offers various
potential benefits, including shorter procedure time, a shorter learning
curve, minimal communication delays and a reduced cognitive load [8]. However,
there are key challenges for effective autonomous RUS. One of the most
important challenge has to do with the hand motions for ultrasound images
acquisition. Such images exhibit considerable inter- and intra-subject
variability and the image quality is highly dependent on the precise position,
orientation and pressure of the ultrasound probe. With incorrect probe
maneuvers, the resulting image presents noise, artifacts, blurred boundaries
and poor visibility, thereby making it unacceptable for diagnosis.
Sonographers rely on visual and haptic feedback, anatomical information, and
diagnostic expertise from prior medical education to rapidly acquire the high-
quality images. Therefore, the RUS must locate the regions with acceptable
diagnostic image quality for inter- and intra-patient procedures in the fewest
exploration steps.

In this paper, we present an autonomous robotic ultrasound system that uses
the domain-expertise in Bayesian Optimization (BO)-based search to scan the
anatomical regions for acquiring diagnostic quality ultrasound images, thereby
eliminating the need to thoroughly scan the entire region. The key
contributions of our work are as follows:

  1. 1.

We proposed a prior in BO, gleaned from the expert’s demonstration of high
image quality probing poses, termed as expert’s prior. BO then estimates the
region’s unknown image quality as a semi-parametric Gaussian process model
with expert’s prior.

  2. 2.

A novel image quality metric is proposed, trained using a dataset of
ultrasound images labelled for diagnostic quality by expert radiologists,
which provides image feedback of the region to the BO.

  3. 3.

We experimentally validated using three urinary bladder phantoms requiring
different probing maneuvers for acquiring high image quality. The results show
that our systems consistently and autonomously acquire high-quality ultrasound
images in all phantoms.

We believe that the use of BO combined with domain expertise to perform
autonomous ultrasound scanning will lead to less reliance on expert
availability and a wider application in remote and underserved populations.

###  I-A Related Work

Autonomous Robotic Ultrasound Systems: In recent years, a range of autonomous
robotic ultrasound systems has been proposed to minimize human intervention.
Earlier works used image features for ultrasound image-based visual servoing
[9, 10, 11]. Later, various systems used pixel-based confidence map methods
[12] and segmentation of structures for optimizing the probe poses and forces
[13, 14, 15, 16]. However, these image feature- and pixel-based approaches are
modality specific, computationally expensive and do not consider the
significance of diagnostic aspects. Hennersperger et al. [17] developed the
autonomous system using the pre-operative MRI scan, however, MRI is quite
expensive to acquire. Ma et al. [18] proposed autonomous lung scanning by
localizing the target region using RGB-D sensor data. However, the system used
only force feedback and did not rely on ultrasound image feedback, thereby
limiting its diagnostic accuracy.

Recently, Li et al. [19, 20] proposed a deep Reinforcement Learning (RL)
framework to control the probe for spinal ultrasound, incorporating image
quality optimization into the reward formulation. However, the success of
these systems is limited to phantoms and patients whose data was included
during training. Moreover, deploying RL in medical systems is quite
challenging, as it requires vast amount of physical interaction with the human
body and poses safety and ethical concerns. In contrast to these systems, the
proposed autonomous ultrasound system narrows down the area to be scanned
using BO, eliminating the need to thoroughly scan the entire region. We
further propose using domain expertise gleaned from the experts in the form of
BO prior and image quality metrics, in order to acquire diagnostic-quality
ultrasound images.

Bayesian Optimization for Medical Robots: Due to the fast optimization
capability, BO has been adopted for safety-critical robotic medical
procedures, such as autonomous robotic palpation [21], semi-autonomous
surgical robot [22], controller tuning of hip exoskeletons [23] and autonomous
robotic ultrasound [24, 25]. Our work is a non-trivial extension to the work
by Goel et al. [25]. They proposed using BO for autonomous ultrasound
utilizing segmentation of the vessel in the ultrasound image as feedback to
the BO for scanning the region with high vessel density. They used hybrid
position-force control to move the robot in $(x,y)$ plane while maintaining
constant force along the $z-$direction to the point of contact. In contrast,
our work suggests two technical improvements to enhance the practicality of
this approach. First, we recommend using a deep learning model that generates
quality scores for ultrasound images as feedback to the BO instead of relying
on a segmented mask of the tissue or structure. The latter approach can be
very time-consuming and labor-intensive for experts as they would need to
annotate anatomical structures’ boundaries, taking into account the ultrasound
image noise and variability due to machine settings, probe pressure, and
patient anatomy. Second, we expand the capabilities of the BO by enabling it
to search for the optimal scanning region along the $(x,y,z)$-axis. Notably,
the $z$-axis is under variable force control to account for varying
physiological conditions [26].

Domain Expertise in BO: BO can utilize the expert’s knowledge in the form of
priors (beliefs) that the expert (practitioner) has on the potential location
of the optimum. Such techniques have been mostly used for hyper-parameter
tuning of image and text datasets [27], open-source machine learning datasets
[28] and robot simulation experiments [29]. A few recent works have utilized
expert’s knowledge in the form of prior for medical robots [30, 31]. Ayvali et
al. [30] propose robotic palpation to detect tissue abnormalities using BO.
They modified the acquisition function of BO, whose value peaks at the user-
provided locations. Zhu et al. [31] proposed an autonomous robotic
auscultation system for locating the optimal sound quality location using BO.
They used visual registration of the patient to locate the anatomical
landmarks for obtaining a prior observation model. Inspired by these works, we
propose BO for autonomous ultrasound leveraging a prior quality map gleaned
from expert’s demonstrations.

Figure 2: Overview of the pipeline for autonomous robotic ultrasound using
online Bayesian optimization (BO), and offline domain expertise to obtain a
prior quality map and to learn image quality assessment metric for providing
feedback to BO.

##  II Methodology

The pipeline of the autonomous robotic ultrasound system is shown in Fig. 2.
In the offline phase, the expert will demonstrate the potential probing poses
to acquire the diagnostic quality images. This demonstrated data would be used
to build a prior quality map, which encodes prior anatomical approximation
about expected image quality. We also built a dataset of urinary bladder
ultrasound images of humans and phantoms with labelled image qualities and
trained a deep learning model for image quality assessment metrics. In the
online phase, we used BO to select the probe poses to find the optimal
ultrasound image quality utilizing both the prior map and quality metric
gleaned from the domain expertise.

###  II-A Bayesian optimization formulation

We use BO to search adaptively for probing poses that yield a high-quality
ultrasound image within a specified anatomical region. Let $A$ be the region
of interest on the human body enclosing the anatomical structure, then the
objective of BO is to solve:

$\max_{\boldsymbol{p}\in A}q(\boldsymbol{\mathcal{I}}(\boldsymbol{p}))$ (1)

where $q(\boldsymbol{\mathcal{I}}(\boldsymbol{p}))$ denotes the quality score
of ultrasound Image $\boldsymbol{\mathcal{I}}$ at probe pose $\boldsymbol{p}$.
The BO will compute the probabilistic estimate of the unknown quality map
$q(\boldsymbol{\mathcal{I}}(\boldsymbol{p}))$ across the human body using the
domain expertise in the form of prior and image quality metric. An acquisition
function is optimized to yield the new probing pose. Once the new observation
is found, the estimate is re-fitted to the data and the process is repeated
till the termination criteria is reached, which is either the maximum
reasonable iteration $N_{max}$ or the estimated quality score threshold
required for adequate diagnosis. The overall algorithm is outlined in
Algorithm 1.

####  II-A1 Expert’s prior

A common estimator used in BO is Gaussian Process (GP) model, which defines an
unknown function $f$ by assigning a probe pose $\boldsymbol{p}$ a random
variable $f(\boldsymbol{p})$, which jointly represent a Gaussian. A GP for
unknown function $f$ is defined by the mean function $\boldsymbol{\mu}(\cdot)$
and covariance or kernel function $\boldsymbol{\kappa}(\cdot,\cdot)$. Given
the function value estimates
$\bar{\boldsymbol{f}}=[{f}(\boldsymbol{p}_{1}),\cdots,{f}(\boldsymbol{p}_{n})]$
at probe poses
$\boldsymbol{\bar{p}}=[\boldsymbol{p}_{1},\cdots,\boldsymbol{p}_{n}]$, GP
regression can predict the function $f$ at new probe pose $\boldsymbol{p}^{*}$
as the Gaussian distribution and is given by:

$\mathcal{P}({f}(\boldsymbol{p}^{*})|\boldsymbol{p}^{*},\bar{\boldsymbol{p}},\bar{\boldsymbol{f}})=\mathcal{N}(\boldsymbol{k}\boldsymbol{K}^{-1}\bar{\boldsymbol{f}},\boldsymbol{\kappa}(\boldsymbol{p}^{*},\boldsymbol{p}^{*})-\boldsymbol{k}\boldsymbol{K}^{-1}\boldsymbol{k}^{T})$
(2)

where,

$\boldsymbol{k}=\begin{bmatrix}\boldsymbol{\kappa}(\boldsymbol{p}_{*},\boldsymbol{p}_{1})&\cdots&\boldsymbol{\kappa}(\boldsymbol{p}_{*},\boldsymbol{p}_{n})\end{bmatrix}$
$\boldsymbol{K}=\begin{bmatrix}\boldsymbol{\kappa}(\boldsymbol{p}_{1},\boldsymbol{p}_{1})&\cdots&\boldsymbol{\kappa}(\boldsymbol{p}_{1},\boldsymbol{p}_{n})\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{\kappa}(\boldsymbol{p}_{n},\boldsymbol{p}_{1})&\cdots&\boldsymbol{\kappa}(\boldsymbol{p}_{n},\boldsymbol{p}_{n})\end{bmatrix}$

We opted to use a combination of two kernel functions, namely the radial basis
function and white noise function, as their combination improved estimations
for structures present in ultrasound images [25]. The formulation of the
kernel is:

${\kappa}(\boldsymbol{p}_{i},\boldsymbol{p}_{j})=\sigma_{r}\exp\bigg{(}\frac{-||\boldsymbol{p}_{i}-\boldsymbol{p}_{j}||^{2}}{2l^{2}}\bigg{)}+\sigma_{w}\boldsymbol{{I}}$
(3)

where $\sigma_{r}$ is the overall variance, $l$ is the length-scale,
$\sigma_{w}$ is the variance of noise and $\boldsymbol{{I}}$ is the identity
matrix. We further denote the set of image qualities as
$\bar{\boldsymbol{q}}=[q_{1},\cdots,{q}_{n}]$.

In GP, we propose using prior knowledge gleaned from expert’s demonstrations
to reduce the explorations and capture the variations of probe poses on the
magnitude of ultrasound image quality corresponding to different human
anatomy. Inspired from work in [31], we formulated the GP as a semi-parametric
GP model, with its prior $\mathcal{E}({\boldsymbol{\theta}})$ modeled as a
Gaussian process with latent parameters $\boldsymbol{\theta}$, representing
the mean $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and covariance function
$\boldsymbol{\kappa}$. The parameters $\boldsymbol{\theta}$ is initially
inferred from observed probe poses and ultrasound image qualities, which the
expert will provide by maneuvering the probe at the potential poses of the
optimum image quality across different subjects. During online BO,
$\boldsymbol{\theta}$ will be inferred using the history of points in
$(\bar{\boldsymbol{p}},\bar{\boldsymbol{q}})$ and prior
$\mathcal{E}({\boldsymbol{\theta}})$ with Maximum A Posteriori (MAP)
estimation, using an L-BFGS solver as:

$\boldsymbol{\theta}^{*}=\arg\max_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta}|\bar{\boldsymbol{p}},\bar{\boldsymbol{q}})\mathcal{E}(\boldsymbol{\theta})$
(4)

where
$\mathcal{L}(\boldsymbol{\theta}|\bar{\boldsymbol{p}},\bar{\boldsymbol{q}})=\prod\mathbb{P}(\boldsymbol{q}_{i}|\boldsymbol{\mu}_{\boldsymbol{\theta}}(\boldsymbol{p}_{i}),\boldsymbol{K})$
is the likelihood function and $\mathbb{P}(.)$ denotes the probability density
function of Gaussian distribution
$\mathcal{N}(\boldsymbol{q}_{i}|\boldsymbol{\mu}_{\boldsymbol{\theta}}(\boldsymbol{p}_{i}),\boldsymbol{K})$.
Since GP models the residual function $f(\boldsymbol{p})$ with respect to the
prior, we subtract the prior from image quality as
$f(\boldsymbol{p}_{i})=q_{i}-{\mu}_{\boldsymbol{\theta}}(\boldsymbol{p})$,
before re-estimating the GP.

Input: Prior $\mathcal{E}{(\boldsymbol{\theta})}$, Region $A$, max. iterations
$N_{max}$;

Initialize $\bar{\mathbf{p}}=\\{\\}$, $\bar{\mathbf{f}}=\\{\\}$,
$\bar{\mathbf{q}}=\\{\\}$;

for  _$i=1,...,N_{max}$ _ do

$\boldsymbol{p}_{i}\leftarrow\arg\max_{\boldsymbol{p}\in
A}EI(\boldsymbol{p})$;

if  _termination criteria reached_ then

stop;

else

Probe at $\boldsymbol{p}_{i}$, compute image quality $q_{i}$;

Set $\bar{\mathbf{p}}\leftarrow\bar{\mathbf{p}}\cup\\{p_{i}\\}$,
$\bar{\mathbf{q}}\leftarrow\bar{\mathbf{q}}\cup\\{q_{i}\\}$;

$\boldsymbol{\theta}\leftarrow\text{argmax}\,\mathcal{L}(\boldsymbol{\theta}|\bar{\mathbf{p}},\bar{\mathbf{q}})\mathcal{E}(\boldsymbol{\theta})$;

Set
$\bar{\mathbf{f}}\leftarrow\bar{\mathbf{f}}\cup\\{q_{i}-{\mu}_{\boldsymbol{\theta}}(\boldsymbol{p})\\}$;

Re-estimate GP;

return Top probe poses with max. image quality;

Algorithm 1 Bayesian Optimization for Ultrasound

####  II-A2 Acquisition Function

In each iteration of BO, the next probe pose to observe the image quality is
determined using an acquisition function. We have used an Expected Improvement
($EI$), which is the most commonly used acquisition function. If the posterior
mean and variance of GP is given by
$\boldsymbol{\mu}_{\bar{\boldsymbol{f}}}(\boldsymbol{x}),\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(x)$,
then $EI$ can be formulated as:

$EI(\boldsymbol{p})=\begin{cases}(\boldsymbol{\mu}_{\bar{\boldsymbol{f}}}(\boldsymbol{p})-f^{+}(\boldsymbol{p})-\xi)\boldsymbol{\Phi}(\boldsymbol{Z})+\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(\boldsymbol{p})\boldsymbol{\phi}(\boldsymbol{Z})&\text{if
}\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(\boldsymbol{p})>0\\\
0&\text{if
}\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(\boldsymbol{p})=0\end{cases}$

(5)

where
$\boldsymbol{Z}=\frac{\boldsymbol{\mu}_{\bar{\boldsymbol{f}}}(\boldsymbol{p})-{f^{+}}(\boldsymbol{p})}{\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(p)}$
if $\boldsymbol{\sigma}_{\bar{\boldsymbol{f}}}^{2}(\boldsymbol{p})>0$ else
$0$; $\boldsymbol{\Phi}$ and $\boldsymbol{\phi}$ are the probability and
cumulative density function of standard normal distribution, respectively and
$f^{+}(\boldsymbol{p})$ is the best observed quality so far. The parameter
$\xi$ in eq. (5) governs the amount of exploration during optimization and a
high $\xi$ value leads to more exploration or less exploitation.

###  II-B Expert’s ultrasound image quality metric

####  II-B1 Dataset

We used two datasets of Urinary Bladder (UB) ultrasound images. One of them is
collected during the in-vivo trials of our in-house developed Telerobotic
Ultrasound System [3, 32] at All India Institute of Medical Sciences (AIIMS),
Delhi, India. The AIIMS ethics committee approved this study under
IEC-855/04.09.2020,RP-16/2020. The other dataset is collected from the UB
phantom. A total of $2016$ real and $2016$ phantom images were collected. The
ground truth quality of the images is an average integer score of labels by
three expert radiologists, each having $15$ years of experience in abdomen
radiology. Each label is an integer score between $1-5$, based on an
internationally prescribed generalized 5-level absolute assessment scale [33,
5]. A score of $1$ means no appearance of the urinary bladder and $5$ means
that the clear depiction of the urinary bladder with distinct boundaries and
acceptable artifacts, depicting a high diagnostic accuracy. A subpar-quality
image ($2$ to $4$) either contains noise or motion artifacts, blurred images,
indistinct boundaries, obscuring the posterior or anterior sections of the
urinary bladder. Later, we normalized the quality score in the range of $0-1$
for standard comparison with other quality estimation methods.

####  II-B2 Feature extraction

Ultrasound image quality assessment requires rich feature extraction for
classifying the images that are highly variable in appearance but differ a lot
in terms of image quality, as shown in Fig. 2. In recent work, Song et al.
[34] proposed a bilinear Convolutional Neural Network (CNN) for fine-grained
classification of breast ultrasound image quality. We propose a technical
enhancement to this work in order to analyze the urinary bladder ultrasound
images, in which the bladder appears at multiple scales/shapes (refer Fig. 2
for sample images) due to the variability among inter- and intra-subject
anatomy, probe poses and forces. Thus, it is also essential to analyze images
at multiple scales. Recently, Basu et al. [35] proposed combining multi-scale
and second-order capability for detecting gall bladder cancer. Taking
inspiration from these works, we proposed a deep CNN-based quality assessment
model. The base network used is Residual Network (ResNet50) [36] and combined
with Multi-scale, Bilinear Pooling classifier, as shown in Fig. 2.

We used the group convolution kernels on equal-width feature volume splits in
place of the $3\times 3$ convolution kernel in the bottleneck layer of
ResNet50. If $\boldsymbol{\mathcal{X}}\in\mathbb{R}^{H\times W\times N}$
represents the feature volume with height $H$, width $W$ and number of
channels $N$, then the operation of Multi-scale block can be represented by
the following equations:

$\displaystyle\boldsymbol{\mathcal{Y}}_{1}$
$\displaystyle=\boldsymbol{\mathcal{X}}_{1}$
$\displaystyle\boldsymbol{\mathcal{Y}}_{3}$
$\displaystyle=\boldsymbol{C}_{2}(\boldsymbol{\mathcal{Y}}_{2}+\boldsymbol{\mathcal{X}}_{3})$
$\displaystyle\boldsymbol{\mathcal{Y}}_{2}$
$\displaystyle=\boldsymbol{C}_{1}(\boldsymbol{\mathcal{X}}_{2})$
$\displaystyle\boldsymbol{\mathcal{Y}}_{4}$
$\displaystyle=\boldsymbol{C}_{3}(\boldsymbol{\mathcal{Y}}_{3}+\boldsymbol{\mathcal{X}}_{4})$

where $\boldsymbol{\mathcal{X}}_{i}\in\mathbb{R}^{H\times W\times N/4}$. Each
split $\boldsymbol{\mathcal{X}}_{i}$ is first concatenated with output of
previous split $\boldsymbol{\mathcal{Y}}_{i-1}$ and then fed to the $3\times
3$ convolutional kernel $\boldsymbol{C}_{i}$ to produce an output
$\boldsymbol{\mathcal{Y}}_{i}$. After passing the image through $16$ multi-
scale blocks, the image feature volume
$\boldsymbol{\mathcal{X}}\in\mathbb{R}^{H\times W\times N}$ is passed through
$1\times 1$ convolution block to reduce the feature volume to
$\boldsymbol{X}\in\mathbb{R}^{H\times W\times N^{{}^{\prime}}}$. Then it is
reshaped to a matrix $\boldsymbol{\mathcal{X}}\in\mathbb{R}^{M\times
N^{{}^{\prime}}}$ where $M=H\times W$. Later, a bilinear pooling is applied
as:

$\displaystyle\boldsymbol{\mathcal{B}}=\frac{1}{N^{{}^{\prime}}}(\boldsymbol{\mathcal{X}}\boldsymbol{\mathcal{X}^{T}})+\epsilon\boldsymbol{I}$
(6)
$\displaystyle\boldsymbol{\mathcal{B}}\xleftarrow[]{}\frac{\boldsymbol{\mathcal{B}}}{||\boldsymbol{\mathcal{B}}||_{2}}\xleftarrow[]{}sign(\boldsymbol{\mathcal{B}})\sqrt{|\boldsymbol{\mathcal{B}}|}$
(7)

where eq. (6) computes the outer product of feature volume, eq. (7) will first
perform the element-wise square-root and then the $l_{2}$ normalisation of the
matrix $\boldsymbol{\mathcal{B}}$. Finally, the flattening of the feature map
is done and then a fully connected layer is utilized to return the ultrasound
image quality score.

###  II-C Robot control

The robot controller will move the probe to the new pose
$\boldsymbol{p}=[x,y,z]$ given by BO, where $(x,y)$ is under position control
and $z$ is under force ($f_{z}$) control. For the safety of phantoms, the
force limits have been set to $20N$ [3]. The orientation of the probe is kept
normal to the point of contact. The hybrid position-force control is used for
controlling the robot. After searching, the robot may execute the top probe
poses with maximum image quality.

##  III Results and Discussions

###  III-A Experimental setup

We conducted the experiments on the laboratory setup of the Robotic Ultrasound
System at Purdue University, USA, consisting of a $7$-DoF Sawyer collaborative
robotic arm (Rethink Robotics, Germany) with Micro Convex MC10-5R10S-3
transducer attached to its end-effector. The US image is captured by the
Telemed Ultrasound machine (Telemed Medical Systems, Italy) and is transferred
to the laptop. The ultrasound was performed on a urinary bladder phantom
(YourDesignMedical, USA). We customized this phantom with the $0.39$ inches
thick (subjected to manual cutting error) rectangular layers of ballistic gel
in order to approximately represent the patient’s body with physiological
differences. Thus, we present our results using three phantoms, termed as
$P0$, $P1$ and $P2$, having $0$, $1$ and $2$ layers, as shown in Fig. 3.

Figure 3: Experimental setup of robotic ultrasound system with three phantoms
of the urinary bladder.

The BO and image quality model have been implemented in Python $3.8$ and
PyTorch $1.11$. ROS has been used to integrate and establish communication
among all components of the setup. For BO Algorithm 1, we used $\xi=0.1$,
$N_{max}=50$, $A\in((0,0.15)m,(0,0.15)m,(8-20)N)$ for $(x,y,f_{z})$. The prior
$\mathcal{E}(\boldsymbol{\theta})$ has been modeled using GP by fitting it to
$10$ potential probing poses and corresponding image qualities.

###  III-B Performance of quality assessment model

We trained the ultrasound image quality assessment model, explained in Section
II-B, using the Categorical Cross Entropy (CCE) as a loss function. We split
the dataset in a $90:10$ ratio as a training and testing dataset. We also used
the transfer learning approach [37] and used the proposed model pre-trained on
ImageNet. The stochastic gradient descent has been used as an optimizer with a
learning rate of $0.005$, momentum of $0.9$ and weight decay of $0.0005$. The
size of the input image to the model is $224\times 224$, batch size is $16$
and the network is trained for $100$ epochs. The results in Table I shows that
the proposed model (ResNet50+MS+BP) achieved an increase in accuracy by
$3.01\%$ on a test set when compared to the ResNet50+BP model proposed in
[34].

TABLE I: Comparison of ultrasound image quality assessment model predictions
on testing set, where accuracy values close to $100$ indicate similarity to
the expert’s quality score.

Image quality score | ResNet50+BP [34] | ResNet50+MS+BP (Proposed)
---|---|---
Precision | Recall | Accuracy | Precision | Recall | Accuracy
$\boldsymbol{1}$ | $92.00$ | $97.87$ | $94.84$ | $93.88$ | $97.87$ | $95.88$
$\boldsymbol{2}$ | $83.33$ | $68.96$ | $75.47$ | $83.87$ | $89.65$ | $86.67$
$\boldsymbol{3}$ | $65.22$ | $75.00$ | $69.77$ | $88.24$ | $75.00$ | $81.08$
$\boldsymbol{4}$ | $93.33$ | $85.71$ | $89.36$ | $91.11$ | $83.67$ | $87.23$
$\boldsymbol{5}$ | $90.48$ | $97.44$ | $93.83$ | $90.48$ | $97.44$ | $93.83$
Average | $\boldsymbol{87.67}$ | $\boldsymbol{87.48}$ | $\boldsymbol{87.34}$ | $\boldsymbol{90.23}$ | $\boldsymbol{90.17}$ | $\boldsymbol{90.05}$

###  III-C Comparing different BO strategies

In order to analyze the effectiveness of the proposed methodology, we have
compared the BO with zero prior to the BO with the proposed expert’s prior. We
illustrated these search strategies using the image feedback having a mean of
the segmented mask of the bladder in the ultrasound image ($q_{S}$) as used in
[25] and having proposed ultrasound image quality metric learned from expert’s
rating ($q_{E}$). For segmentation, we used a U-net-based segmentation model
proposed in [38]. Further, each of the feedback strategies has been compared
with different search spaces, first considering the probe motion along $x$ and
$y$-axis of the phantom, second along $x$, $y$ and $z-$axis of the phantom,
where $z-$axis is under the force control $(f_{z})$. The estimated quality
maps obtained using these strategies for $P0$ are shown in Fig. 4, where red
region shows the high-quality region and blue region shows the low-quality
region. The black dots over the map represents the queried probe positions
over the phantom during the optimization. The first column in Fig. 4 shows the
quality map obtained using the uniform movement of the probe over the phantom,
which has been considered as the approximate ground truth quality map. For
both the quality types, the ground truth has been obtained using the
approximate desired force ($f_{d}$) of $14N$, $16N$ and $18N$ for $P0$, $P1$
and $P2$, respectively, which gives the best image quality in these phantoms.
We present results for $3$ cases to illustrate the effect of searching with
appropriate force in these phantoms: (i) $f_{z}<f_{d}$: $fz$ is constant but
equal to $f_{d}-4$, (ii) $f_{z}=f_{d}$ and (iii) when $f_{z}$ is variable.

We compared the quality maps of these strategies by doing quantitative
analysis using three metrics: (i) Sum of quality difference of top $n$ points,
(ii) Top quality, and (iii) Zero Normalized Cross Correlation (ZNCC), as shown
in Table II. The numbers in the table represent the average value of the
matrices for the $3$ tests on each phantom. These metrics have been computed
with respect to the approximate ground truth for the phantom. The sum of the
difference between the top $n-$points compares the quality of images acquired
from the top-n highest quality values, top quality compares the highest value
of image quality score and ZNCC evaluates the overall similarity of the
acquired quality map during the search. The value of quality differences close
to $0$, and top quality and ZNCC value close to $1$ indicates a better
estimation of the quality map. The quality maps in Fig. 4 with less scattered
probe points (less exploration) and more points in the high-quality region
(red) represent a better search strategy.

Figure 4: The estimated ultrasound image quality map of urinary bladder
phantom $P0$ using different BO strategies. Black dots are positions where
probe evaluated the quality. The corresponding ultrasound images are available
in the attached media. TABLE II: Quantitative comparison of different BO
strategies for three different urinary bladder phantoms $P0,P1$ and $P2$

Image quality estimation method | Variables | BO with zero prior | BO with expert’s prior
---|---|---|---
Sum of quality difference of $n$ points | Top | ZNCC | Sum of quality difference of $n$ points | Top | ZNCC
$n=1$ | $n=5$ | $n=10$ | $n=20$ | Quality | $n=50$ | $n=1$ | $n=5$ | $n=10$ | $n=20$ | Quality | $n=50$
Segmentation | $x,y,f_{z}<f_{d}$ | $0.382$ | $0.944$ | $1.430$ | $1.931$ | $0.684$ | $0.689$ | $0.291$ | $0.692$ | $0.941$ | $1.205$ | $0.782$ | $0.782$
($q_{S}$) | $x,y,f_{z}=f_{d}$ | $0.132$ | $0.609$ | $0.963$ | $1.651$ | $0.911$ | $0.811$ | $0.103$ | $0.531$ | $0.785$ | $1.308$ | $0.911$ | $0.920$
[38] | $x,y,f_{z}$ | $0.396$ | $1.016$ | $1.919$ | $3.336$ | $\mathbf{0.711}$ | $\mathbf{0.733}$ | $0.404$ | $0.991$ | $1.379$ | $2.158$ | $\mathbf{0.799}$ | $\mathbf{0.801}$
Expert’s image | $x,y,f_{z}<f_{d}$ | $0.280$ | $0.370$ | $0.600$ | $1.030$ | $0.690$ | $0.717$ | $0.130$ | $0.290$ | $0.570$ | $1.010$ | $0.750$ | $0.817$
quality metric | $x,y,f_{z}=f_{d}$ | $0.120$ | $0.240$ | $0.390$ | $0.710$ | $0.950$ | $0.876$ | $0.050$ | $0.090$ | $0.180$ | $0.970$ | $0.980$ | $0.959$
($q_{E}$) | $x,y,f_{z}$ | $0.130$ | $0.270$ | $0.820$ | $1.320$ | $\mathbf{0.823}$ | $\mathbf{0.821}$ | $0.040$ | $0.280$ | $0.760$ | $1.600$ | $\mathbf{0.910}$ | $\mathbf{0.889}$

From the result in Fig. 4 and Table II, it has been found that the BO using
the segmented image as quality score in $(x,y)$ space with $f_{z}\leq f_{d}$
have resulted in being too exploratory (low ZNCC) with a lot of points spread
over the low-quality region of the phantom. However, the quality maps obtained
using the expert’s quality metric of the image have fewer explorations, with
most of the probe positions in the high-quality region of the phantom. Due to
noise and shadows in the ultrasound image, the segmentation results are prone
to errors, resulting in a large number of probe evaluations in low-quality
regions, whereas expert’s image quality score, which is based on the holistic
assessment of the image, pinpoints the focus on anatomical structures rather
than getting distracted by noise. The search strategies using $fz<f_{d}$ could
not find the high-quality region and instead converged to the local maxima
rather than the global maxima. However, with $f_{z}=f_{d}$, the high-quality
regions have been acquired. When the quality region is searched using $f_{z}$
as a variable in BO with zero prior, the quality maps and top quality score
show that the high-quality regions can be located with a varying force too,
which is essential for in-human ultrasound procedures. However, the search is
quite exploratory, reporting low ZNCC values of $0.733$ and $0.821$ for
quality $q_{S}$ and $q_{E}$, respectively. When the expert’s prior is used,
all BO strategies have significantly improved, including the search space with
three variables ($x,y,f_{z}$). The exploration steps of BO usually increase as
the search space dimension increases. However, BO with expert’s prior reported
a top quality of $0.910$ with a ZNCC score of $0.889$, which is $9.6\%$ and
$7.6\%$ more than the BO with zero prior.

###  III-D Validating the convergence of probe positions and forces

Since our study involves phantom experiments, the approximate probe positions
and forces that yield the best-quality images are known. The search strategy
should converge to these approximate probe poses and forces to acquire high-
quality images. The proposed strategy has reached the desired probe position
with an average mean value accuracy of $98.73\%$ across all phantoms. To
emphasize the convergence of force, we compared the probe forces explored by
different BO search strategies, as shown in Fig. 5.

Figure 5: Force $f_{z}$ profiles with different BO search strategies

The proposed formulation of BO using the expert’s prior and image quality
metric has resulted in the mean value accuracy of $99.28\%$, $98.25\%$, and
$96.11\%$ for $P0$, $P1$, and $P2$, respectively. Comparatively, the other BO
search strategies using zero-prior and segmentation-based quality maps
($q_{S}$) have shown significant errors in mean values and greater standard
deviation due to the noise in image feedback and the inability to adapt to the
profile of the scanning region.

##  IV Conclusion

We proposed an autonomous Robotic Ultrasound System (RUS) to perform the
ultrasound as per clinical protocols. We used Bayesian Optimization (BO) to
search for high-quality regions leveraging the domain expertise in the form of
a prior quality map and ultrasound image quality. The prior map has been
gleaned using expert’s demonstration of the potential high-quality probing
maneuvers. A novel image quality metric has been learned from the expert-
labelled dataset of ultrasound images. Three phantom experiments validated
that incorporating domain expertise into BO effectively improves the system
performance, resulting in acquiring diagnostic quality ultrasound images while
adapting to desired probing maneuvers. Since phantom results are promising, we
would like to validate its capability for in-vivo study using our RUS in India
[3], which is our future work. We would also expand the search space in BO
from $[x,y,f_{z}]$ to include $[roll,pitch,yaw]$ in order to orient the probe
for scanning patients with complex physiological conditions.

## References

  * [1] A. Rykkje, J. F. Carlsen, and M. B. Nielsen, “Hand-held ultrasound devices compared with high-end ultrasound systems: a systematic review,” Diagnostics, vol. 9, no. 2, p. 61, 2019.
  * [2] S. J. Adams, B. Burbridge, H. Obaid, G. Stoneham, P. Babyn, and I. Mendez, “Telerobotic sonography for remote diagnostic imaging: narrative review of current developments and clinical applications,” Journal of Ultrasound in Medicine, vol. 40, no. 7, pp. 1287–1306, 2021.
  * [3] D. Raina, H. Singh, S. K. Saha, C. Arora, A. Agarwal, S. Chandrashekhara, K. Rangarajan, and S. Nandi, “Comprehensive telerobotic ultrasound system for abdominal imaging: Development and in-vivo feasibility study,” in 2021 International Symposium on Medical Robotics (ISMR), pp. 1–7, IEEE, 2021\.
  * [4] J. Wang, C. Peng, Y. Zhao, R. Ye, J. Hong, H. Huang, and L. Chen, “Application of a robotic tele-echography system for covid-19 pneumonia,” Journal of Ultrasound in Medicine, vol. 40, no. 2, pp. 385–390, 2021.
  * [5] S. Duan, L. Liu, Y. Chen, L. Yang, Y. Zhang, S. Wang, L. Hao, and L. Zhang, “A 5g-powered robot-assisted teleultrasound diagnostic system in an intensive care unit,” Critical Care, vol. 25, no. 1, pp. 1–9, 2021.
  * [6] S. J. Adams, B. E. Burbridge, A. Badea, L. Langford, V. Vergara, R. Bryce, L. Bustamante, I. M. Mendez, and P. S. Babyn, “Initial experience using a telerobotic ultrasound system for adult abdominal sonography,” Canadian Association of Radiologists’ Journal, vol. 68, no. 3, pp. 308–314, 2017.
  * [7] E. Rojas-Muñoz and J. P. Wachs, “Assessing task understanding in remote ultrasound diagnosis via gesture analysis,” Pattern Analysis and Applications, vol. 24, no. 4, pp. 1489–1500, 2021.
  * [8] K. Li, Y. Xu, and M. Q.-H. Meng, “An overview of systems and techniques for autonomous robotic ultrasound acquisitions,” IEEE Transactions on Medical Robotics and Bionics, vol. 3, no. 2, pp. 510–524, 2021.
  * [9] R. Mebarki, A. Krupa, and F. Chaumette, “2-d ultrasound probe complete guidance by visual servoing using image moments,” IEEE Transactions on Robotics, vol. 26, no. 2, pp. 296–306, 2010.
  * [10] R. Nakadate, J. Solis, A. Takanishi, E. Minagawa, M. Sugawara, and K. Niki, “Out-of-plane visual servoing method for tracking the carotid artery with a robot-assisted ultrasound diagnostic system,” in 2011 IEEE International Conference on Robotics and Automation, pp. 5267–5272, IEEE, 2011\.
  * [11] C. Nadeau, A. Krupa, J. Petr, and C. Barillot, “Moments-based ultrasound visual servoing: From a mono-to multiplane approach,” IEEE Transactions on Robotics, vol. 32, no. 6, pp. 1558–1564, 2016.
  * [12] A. Karamalis, W. Wein, T. Klein, and N. Navab, “Ultrasound confidence maps using random walks,” Medical image analysis, vol. 16, no. 6, pp. 1101–1112, 2012.
  * [13] P. Chatelain, A. Krupa, and N. Navab, “Confidence-driven control of an ultrasound probe,” IEEE Transactions on Robotics, vol. 33, no. 6, pp. 1410–1424, 2017.
  * [14] Z. Jiang, M. Grimm, M. Zhou, J. Esteban, W. Simson, G. Zahnd, and N. Navab, “Automatic normal positioning of robotic ultrasound probe based only on confidence map optimization and force measurement,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 1342–1349, 2020.
  * [15] Z. Jiang, Z. Li, M. Grimm, M. Zhou, M. Esposito, W. Wein, W. Stechele, T. Wendler, and N. Navab, “Autonomous robotic screening of tubular structures based only on real-time ultrasound imaging feedback,” IEEE Transactions on Industrial Electronics, vol. 69, no. 7, pp. 7064–7075, 2021\.
  * [16] Z. Jiang, Y. Gao, L. Xie, and N. Navab, “Towards autonomous atlas-based ultrasound acquisitions in presence of articulated motion,” IEEE Robotics and Automation Letters, 2022.
  * [17] C. Hennersperger, B. Fuerst, S. Virga, O. Zettinig, B. Frisch, T. Neff, and N. Navab, “Towards mri-based autonomous robotic us acquisitions: a first feasibility study,” IEEE transactions on medical imaging, vol. 36, no. 2, pp. 538–548, 2016.
  * [18] X. Ma, Z. Zhang, and H. K. Zhang, “Autonomous scanning target localization for robotic lung ultrasound imaging,” in 2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 9467–9474, IEEE, 2021\.
  * [19] K. Li, J. Wang, Y. Xu, H. Qin, D. Liu, L. Liu, and M. Q.-H. Meng, “Autonomous navigation of an ultrasound probe towards standard scan planes with deep reinforcement learning,” in 2021 IEEE International Conference on Robotics and Automation (ICRA), pp. 8302–8308, IEEE, 2021.
  * [20] K. Li, Y. Xu, J. Wang, D. Ni, L. Liu, and M. Q.-H. Meng, “Image-guided navigation of a robotic ultrasound probe for autonomous spinal sonography using a shadow-aware dual-agent framework,” IEEE Transactions on Medical Robotics and Bionics, vol. 4, no. 1, pp. 130–144, 2021.
  * [21] Y. Yan and J. Pan, “Fast localization and segmentation of tissue abnormalities by autonomous robotic palpation,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 1707–1714, 2021.
  * [22] J. Chen, D. Zhang, A. Munawar, R. Zhu, B. Lo, G. S. Fischer, and G.-Z. Yang, “Supervised semi-autonomous control for surgical robot based on banoian optimization,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2943–2949, IEEE, 2020.
  * [23] Y. Ding, M. Kim, S. Kuindersma, and C. J. Walsh, “Human-in-the-loop optimization of hip assistance with a soft exosuit during walking,” Science robotics, vol. 3, no. 15, p. eaar5438, 2018.
  * [24] Y. Huang, W. Xiao, C. Wang, H. Liu, R. Huang, and Z. Sun, “Towards fully autonomous ultrasound scanning robot with imitation learning based on clinical protocols,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 3671–3678, 2021.
  * [25] R. Goel, F. Abhimanyu, K. Patel, J. Galeotti, and H. Choset, “Autonomous ultrasound scanning using bayesian optimization and hybrid force control,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 8396–8402, IEEE, 2022.
  * [26] M. Akbari, J. Carriere, T. Meyer, R. Sloboda, S. Husain, N. Usmani, and M. Tavakoli, “Robotic ultrasound scanning with real-time image-based force adjustment: quick response for enabling physical distancing during the covid-19 pandemic,” Frontiers in Robotics and AI, vol. 8, p. 645424, 2021\.
  * [27] Z. Wang, G. E. Dahl, K. Swersky, C. Lee, Z. Mariet, Z. Nado, J. Gilmer, J. Snoek, and Z. Ghahramani, “Pre-training helps bayesian optimization too,” arXiv preprint arXiv:2207.03084, 2022.
  * [28] C. Hvarfner, D. Stoll, A. Souza, M. Lindauer, F. Hutter, and L. Nardi, “pibo: Augmenting acquisition functions with user beliefs for bayesian optimization,” arXiv preprint arXiv:2204.11051, 2022.
  * [29] Z. Wang, B. Kim, and L. P. Kaelbling, “Regret bounds for meta bayesian optimization with an unknown gaussian process prior,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [30] E. Ayvali, A. Ansari, L. Wang, N. Simaan, and H. Choset, “Utility-guided palpation for locating tissue abnormalities,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 864–871, 2017.
  * [31] Y. Zhu, A. Smith, and K. Hauser, “Automated heart and lung auscultation in robotic physical examinations,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 4204–4211, 2022.
  * [32] S. H. Chandrashekhara, K. Rangarajan, A. Agrawal, S. Thulkar, S. Gamanagatti, D. Raina, S. K. Saha, and C. Arora, “Robotic ultrasound: An initial feasibility study,” World Journal of Methodology, vol. 12, no. 4, pp. 274–284, 2022.
  * [33] P. Cantin and K. Knapp, “Diagnostic image quality in gynaecological ultrasound: Who should measure it, what should we measure and how?,” Ultrasound, vol. 22, no. 1, pp. 44–51, 2014.
  * [34] Y. Song, Z. Zhong, B. Zhao, P. Zhang, Q. Wang, Z. Wang, L. Yao, F. Lv, and Y. Hu, “Medical ultrasound image quality assessment for autonomous robotic screening,” IEEE Robotics and Automation Letters, vol. 7, no. 3, pp. 6290–6296, 2022.
  * [35] S. Basu, M. Gupta, P. Rana, P. Gupta, and C. Arora, “Surpassing the human accuracy: Detecting gallbladder cancer from usg images with curriculum learning,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pp. 20886–20896, 2022.
  * [36] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770–778, 2016.
  * [37] P. M. Cheng and H. S. Malhi, “Transfer learning with convolutional neural networks for classification of abdominal ultrasound images,” Journal of digital imaging, vol. 30, no. 2, pp. 234–243, 2017.
  * [38] D. Raina, K. Verma, S. Chandrashekhara, and S. K. Saha, “Slim u-net: Efficient anatomical feature preserving u-net architecture for ultrasound image segmentation,” in 2022 International Conference on Biomedical and Bioinformatics Engineering (ICBBE), pp. 1–8, ACM, 2022.