\section{Introduction}
The Kondo effect~\cite{Kondo1964}~ is responsible for the low temperature resistivity upturn in metals with dilute magnetic impurities~\cite{DeHaas1934}~, and is a paradigmatic problem in Condensed Matter Physics.
In the Kondo effect, the coupling of a single spin to a metallic environment causes the formation of a magnetic screening cloud.
It embodies the complexity arising when a single particle interacts with a many-body environment, and was studied extensively using electronic transport through GaAs quantum dots (QDs)~\cite{Goldhaber-Gordon1998}~, and later carbon nanotubes~\cite{Nygard2000}~ and single molecules~\cite{Park2002}~. In these systems, the effect manifests itself as an enhancement of the electronic conductance due to correlations between the source and drain mediated by an unpaired spin localized within the QD.

A natural extension of the Kondo problem involves coupling a QD to one or two superconducting leads.
The first experimental study where a QD in the Kondo regime was coupled to superconducting (S) source and drain electrodes was reported by Buitelaar et al. ~\cite{Buitelaar2002}~. In this `S-QD-S' geometry, the electron number in the impurity level was tuned to an odd value via electronic gating. The SU(2) Kondo signature – zero-energy state - was observed whenever the Kondo temperature $T_K$ became the dominant energy scale in the system.
Such an observation falls in line with the theory of many-body singlet state for the Kondo ground state, which can create a Kondo signature whenever the binding energy of the Kondo singlet exceeds the superconducting energy gap $\Delta$~\cite{Clerk2000}~.

Later studies have shown that even when the zero-bias Kondo feature is  suppressed, an applied source-drain bias $V_{SD} = \Delta$ reveals the formation of a Kondo resonance within the QD - the strength of which appears to follow a scaling law defined by the competition between $\Delta$ and $k_{B}T_{K}$~\cite{Buizert2007,Clerk2000}~, where $k_B$ is the Boltzmann constant.
This competition underlies a rich set of physical phenomena involving the Andreev reflection at the QD-S interface and the Kondo singlet formation at the QD-N (normal) interface ~\cite{Clerk2000,Yeyati2003,Sand-Jespersen2007}~. Theoretical studies have suggested that adopting a hybrid N-QD-S geometry would allow a more rigorous approach in studying these dynamics ~\cite{Clerk2000,Cuevas2001,Li2016,Zitko2015}~.

Experimental realizations of N-QD-S devices typically rely on clean semiconductor nanowires (GaAs, InAs)~\cite{Jespersen2006}~, coupled to superconductor materials that could be evaporated as electrodes (Pb, Nb, Al, MoRe).
The advent of van der Waals (vdW) materials now enables the study of a whole library of materials with diverse ground states ~\cite{Ajayan2016}~. Specifically, the Transition Metal Dichalcogenides (TMD) family contains candidate 2D superconductors (NbSe$_2$\xspace, NbS$_2$\xspace, PdTe$_2$\xspace) that exhibit superconductivity from the bulk regime to the monolayer limit ~\cite{Ugeda2016}~.
Insulators such as MoS$_2$\xspace, \text{WSe$_2$}, and hexagonal Boron Nitride (hBN) have been successfully implemented as barriers in tunnel junctions ~\cite{Dvir2018a, Dvir2018b, khestanova2018}~ and have been shown to host atomic defects ~\cite{Zhou2013, Dvir2019, Devidas2021}~ which could be used as QDs ~\cite{Chandni2015, Greenaway2018, Keren2020}~.

In this work we demonstrate that a Kondo effect can also be realized by coupling a QD to a layered superconductor. To place a QD in close proximity to a layered superconductor, we make use of naturally-occurring defects in a TMD semiconductor tunnel barrier placed on top of the superconductor NbSe$_2$\xspace. As we have shown in previous studies, such defect-bound dots may couple strongly to the SC, giving rise to Andreev bound state (ABS) sub-gap conductance features~\cite{Dvir2019}~. Conversely, when the QD is weakly coupled, it can serve as a sensitive spectral~\cite{Devidas2021}~ or compressibility~\cite{Keren2020}~ probe.
In the present case, we find a tell-tale zero-bias conductance peak which sets in at finite in-plane and out-of-plane magnetic field. The peak splits at higher magnetic fields, consistent with a Kondo feature.
In NbSe$_2$\xspace, tunneling measurements resolve two gaps~\cite{Noat2015,Dvir2018a}~. Here we find that the Kondo feature is correlated with a suppression in the spectral signature of the lower of these two energy gaps, associated with the Se band.
The correlation of the peak value with the suppression of the second gap tunneling feature suggests that the Kondo effect is mediated by carriers belonging to the Se-derived band, which turns normal at a low magnetic field, while the larger Nb-derived gap remains stable. Our results suggest that defect-bound QDs can be used to study the Kondo effect in unconventional, van der Waals, layered superconductors.

\section{Experimental Details}
The schematic of the MoS$_2$\xspace-graphene-NbSe$_2$\xspace vdW heterostructure is shown in Figure \ref{fig:fig1}(a) and the optical image of the measured device reported in this work is shown in Figure \ref{fig:fig1}(b).
MoS$_2$\xspace, graphene and NbSe$_2$\xspace are exfoliated independently on Si/SiO$_2$ (285 nm oxide) substrates. The desired flakes of 3-4 layers MoS$_2$\xspace (tunnel barrier), monolayer graphene and bulk NbSe$_2$\xspace are chosen by optical contrast. The MoS$_2$\xspace barrier flake is first picked up using the polycarbonate (PC) technique. Graphene is subsequently picked up using the van der Waals interaction between MoS$_2$\xspace and graphene. The picked up heterostructure is then transferred on to a 25 nm thick bulk NbSe$_2$\xspace flake. The process is carried out in a glovebox under an Argon atmosphere. Standard e-beam lithographic techniques are used to pattern tunnel electrodes and ohmic contacts on MoS$_2$\xspace and NbSe$_2$\xspace respectively. Ti/Au electrodes are evaporated using an e-beam evaporator with an additional Argon ion milling step prior to the ohmic contact deposition step, so as to obtain better contact on the NbSe$_2$\xspace flake. No contacts are made to the graphene flake.

\section{Results}
\begin{figure}
    \centering
    \includegraphics[width = 400 pt]{figure1_v7.png}
    \caption{(a) Schematic of the device measured in the current work; (b) Optical image of the measured device. NbSe$_2$\xspace, graphene and MoS$_2$\xspace flakes are outlined and labelled in different colours for easier identification. S and D indicate the source and drain electrodes respectively. (c) Left $y$-axis - Differential conductance ($dI/dV$) spectrum of the NbSe$_2$\xspace-graphene-MoS$_2$\xspace tunnel junction depicted in (b), Right $y$-axis - Derivative of the differential conductance ($d^2I/dV^2$) spectrum measured at 28 mK. The black dashed lines track the two NbSe$_2$\xspace gaps present in the spectrum at 0.6 mV and 1.25 mV; (d) and (e) show the evolution of the tunneling spectrum as a function of parallel ($B_{||}$) and perpendicular magnetic field ($B_{\perp}$) respectively, for low magnetic field values. }
   \label{fig:fig1}
\end{figure}

The tunneling differential conductance ($dI/dV$) of the MoS$_2$\xspace-graphene-NbSe$_2$\xspace\ stack at 28 mK and B = 0 T, is shown on the left $y$-axis of Figure 1(c). The spectrum exhibits well-defined quasiparticle peaks and a hard gap - manifest in a ratio of over 100 between the conductance outside the gap ($G_N$) and zero-bias conductance ($G_0$). Such a hard gap, seen in our earlier studies on NbSe$_2$\xspace tunneling devices~\cite{Dvir2018a,Dvir2018b}~, attests to the quality of the tunnel junction, which suppresses two-particle tunneling contributions.
The $2^{nd}$ derivative of the tunneling current ($d^2I/dV^2$), plotted in the right $y$-axis of Figure 1(c) shows a clear separation between two distinct features marked by black dashed lines. These two features correspond to the higher energy gap, $\Delta_1 = 1.25$ mV, and the lower energy gap $\Delta_2 = 0.6$ mV. These were identified in the past as related to the two superconducting bands in NbSe$_2$\xspace – the 2D-like niobium $d$-orbitals bands around the $\Gamma, K$ points, associated with $\Delta_1$, and selenium-derived $p$-orbitals which form a small 3D Fermi surface around the $\Gamma$-point, associated with $\Delta_2$ ~\cite{Dvir2018a,Dvir2018b,Kiss2007,Noat2015}~.

Interestingly, the tunneling spectrum does not reveal any signature which can be associated with the graphene layer.
This can be understood considering that only for a small range of twist angles graphene's Fermi pockets
can have some overlap with the Fermi pockets of a TMD like NbSe$_2$\xspace~\cite{Gani_NbSe2_Gr_2019}~.
As a result, the graphene spectrum may be decoupled from the NbSe$_2$\xspace.

Figure 1(d) shows the evolution of the tunneling spectra as a function of in-plane magnetic field ($B_\parallel$) applied parallel to the sample plane. The observed spectra evolve in a way which is not characteristic of NbSe$_2$\xspace devices measured in the past ~\cite{Dvir2018a}~: First, we find that the superconducting gap is filled even by the application of a very low field. Second, a zero-bias conductance peak (ZBCP) develops at $B_\parallel = 0.275$ T. In what follows, we discuss possible mechanisms for these features.

Zero-bias spectral features are rather common in proximity superconducting devices, and may have a number of possible origins.
They could be associated with Andreev bound states (ABS) residing on the surface of a $d$-wave superconductor ~\cite{Millo2018}~, or with constructive superposition of bound states formed by reflectionless tunneling at diffusive N-S junctions ~\cite{VanWees1992,Marmorkos1993}~.

Zero-bias states may appear when a QD, proximity-coupled to  a superconductor, forms a local Andreev bound state which undergoes a singlet-doublet transition~\cite{Lee2014,Scherubl2020}~.
Conversely, when the QD is weakly coupled to both SC ($\Gamma_{SC}$) and normal ($\Gamma_N$) leads, i.e.
$\Gamma_{SC,N} << \Delta$,
QD-SC transport is dominated by single electron resonant tunneling.
In this regime, the spectrum exhibits sharp conductance peaks, but these can appear both above the gap and below it,
when the density of states (DOS) is not zero.
In a recent publication~\cite{Devidas2021}~ we have reported the use of such a weakly-coupled QD as a sensitive spectrometer.

The zero-bias conductance peak also appears as a response to out-of-plane magnetic field $B_{\perp}$ as seen in Figure 1(e). To further study the nature of this peak, we track it's evolution as magnetic fields are applied in directions parallel and perpendicular to the sample plane.
Figure 2(a) represents a color scale map of the tunneling spectra vs. $B_\parallel$ extending to 9 T. Following the evolution of the zero-energy state that switches on at 0.275 T (Figure 1(d)), the peak persists at the same energy until a magnetic field induced splitting of the state into two distinct features is observed at 1 T.
The splitting is symmetrical in $V_{SD}$ with respect to zero.

Although such linear dispersions are regularly observed in ABS, here we argue that these features are unlikely to be related with ABS, on two grounds. First, ABS features in NbSe$_2$\xspace usually appear at finite bias, and may converge to zero at finite field. They are observed starting at zero magnetic field and do not require a finite field to be visible. Second, ABS features in NbSe$_2$\xspace tunnel devices are observed only in very thin flakes. In any flake thicker than a few layers, the in-plane magnetic field introduces a sub-gap tunneling signal which obscures them ~\cite{Dvir2019}~.

We are thus led to propose that the observed feature - the zero-bias conductance peak at a finite field, is associated with a Kondo origin.

The split states are separated by 2$E_Z$ ~\cite{Goldhaber-Gordon1998,Pustilnik2004}~, $E_Z=\pm g\mu_B B$ being the Zeeman energy, with a Land\'e $g$-factor (1.67$\pm$0.04).
A similar trend is observed in the evolution of the spectrum in $B_\perp$ with a field induced splitting appearing at 0.6 T, shown in Figure 2(b). The value of Land\'e $g$-factor obtained from out-of-plane data is 1.94$\pm$0.04. We note that similar values have been seen in the past for atomic defect QDs~\cite{Dvir2019,Devidas2021}~.

\begin{figure}
    \centering
    \includegraphics[width = 420 pt]{figure2_v6.png}
    \caption{Color map of the $dI/dV$ spectrum in (a) in-plane magnetic field ($B_{||}$) and (b) out-of-plane magnetic field ($B_\perp$). The black arrows in (a) and (b) indicate the Zeeman split features.}
    \label{fig:fig2}
\end{figure}

The Kondo effect sets in below a Kondo temperature $T_K$ which can be obtained through temperature dependence or through the profile of the zero-bias conductance peak.
Temperature dependence measurements of the zero-bias conductance peak ($B_\perp$ = 25 mT) are carried out from 28 mK to 419 mK, beyond which the peak is no longer detectable above the background. A selection of these plots is shown in Figure 3(a). The maximal zero-bias conductance at 28 mK is 11.1 $\mu S$ = 0.14 $e^2/h$. Figure 3(b) shows a semi-log plot of zero-bias conductance peak heights as a function of temperature, where peak height values are obtained after suitable background subtraction.
The data is fit to an empirical equation \eqref{eq1} derived from the numerical renormalization group (NRG) theory for the Kondo ground state~\cite{Buitelaar2002,Goldhaber-Gordon1998}~.

\begin{equation}
\label{eq1}
\centering
    G(T)=\frac{G_{max}}{\bigg[1+\bigg (2^{\nicefrac{1}{s}}-1\bigg)\bigg(\frac{T}{T_K} \bigg)^2 \bigg]^s}
\end{equation}

Where $G_{max}$ is the maximum conductance at the lowest temperature measured, $T_K$ is the Kondo temperature and $s$ is a dimensionless value related to the spin-state of the electron in the QD. The best fit to the data (solid red line) yields $G_{max}$ = 0.0695$\pm$0.003 $e^2/h$; $T_K$ = 153$\pm$6 mK (13.2 $\mu$eV) and $s$ = 1.4$\pm$0.3.

Clearly, the exponent $s$ deviates from the value $s$ = 0.22, expected for a spin 1/2 system at the Kondo regime~\cite{Goldhaber-Gordon1998,Buitelaar2002}~. Values removed from 0.22 could indicate that the QD is in the mixed-valence regime, which is seen when it's energy is close to resonance, or, alternatively, that the spin is not 1/2.

\begin{figure}[ht]
    \centering
    \includegraphics[width = 400 pt]{figure3_v12.png}
    \caption{Evolution of the sub-gap conductance ($dI/dV$) spectra as a function of temperature in the QD-NbSe$_2$\xspace junction (a) Scans at select temperatures; (b) Kondo peak heights as a function of temperature, after suitable background subtraction, fitted (red solid line) to the NRG derived empirical equation  ~\cite{Buitelaar2002,Goldhaber-Gordon1998}~; (c) The zero-bias conductance peak at 28 mK, fit to the Fano function ~\cite{Fano1961}~ to extract the value of HWHM which defines $T_K$ ~\cite{Deacon2010b}~; (d) The zero-bias  conductance peak at 28 mK, fit to the equation derived by Kretinin et al. ~\cite{Kretinin2012}~ to different values of $T_K$ and $s$.}
    \label{fig:fig3}
\end{figure}

It is also possible to evaluate $T_K$ from the conductance feature line-shape.
In a N-QD-S system, the half-width at half maximum (HWHM) of the peak at the lowest temperature is directly related to $T_K$ ~\cite{Deacon2010b}~. Since the Kondo state involves an interaction between a discrete state (QD) and continuum (N lead), we use the Fano function ~\cite{Fano1961}~ to extract the HWHM (Figure 3(c)). The function is defined as

\begin{equation}
\label{eq_Fano}
\centering
    G(V_{SD})= A \frac{(\varepsilon+q)^2}{(1+\varepsilon^2 )}+B; \varepsilon= \frac{(V_{SD}-\varepsilon_0 )}{\Gamma_d}
\end{equation}

Where $\varepsilon_0$ is the resonance energy (0 eV in the present data), A and B are constants, $\Gamma_d$ is the HWHM, and $q$ is a phenomenological dimensionless shape parameter. The HWHM obtained is $\Gamma_d$ = 25.8$\pm$6 $\mu$eV which corresponds to $T_K$ = 300$\pm$70 mK.
The other parameters obtained from the fit are A = -5.5$\pm$0.7  $\mu$S, B = 0.11$\pm$0.06 $\mu$S, q = 0.10$\pm$0.05.

The ambiguity in the extracted Kondo temperature has been observed earlier on N-QD-N devices too ~\cite{Kretinin2012,VanderWiel2000}~. Kretinin et al. ~\cite{Kretinin2012}~ improved on the earlier NRG formula by suggesting that analysing the zero energy state at the lowest temperature as a function of a normalised energy scale $\nu\equiv eV_{SD}/(k_BT_K )$ is a more reliable method to extract the correct Kondo temperature of the system. It was based on the argument that temperature dependence data might have additional features arising from non-Kondo origin incorporated in the zero-energy state as a result of which the extracted $T_K$ would not define the Kondo state alone. The argument might hold for the anomaly in the value of $s$ parameter we obtained using the empirical formula fit.

Using the equation derived by Kretinin et al.:
\begin{equation}
\label{eq_Kretinin}
\centering
    G(T=0,\nu)= \frac{G_{max}} {\bigg[1+\bigg(2^{\nicefrac{1}{s}}-1\bigg)\frac{\nu^2}{\pi} \bigg]^{s}}
\end{equation}
and assuming our base temperature is close to the condition T = 0, we plot our 28 mK data (Figure 3(d)) along with two curves generated using Equation~\eqref{eq_Kretinin}. We notice that the curve generated by the $T_K$ and $s$ values obtained from the temperature dependence fit (blue dash dot) doesn't replicate the behaviour of the measured data. However, for the values $T_K$ = 270 mK (23.26 $\mu$eV) and $s$ = 0.22 (red dash dot), the curve closely follows the measured data.

\section{Discussion}
The results presented so far suggest that the zero-bias conductance peak we observe is indeed associated with a Kondo feature. Yet the larger gap $\Delta_1$ = 1.25 mV is far greater than $k_BT_K$, regardless of the method we choose to evaluate $T_K$.
Although theoretical models do suggest that such a co-existence is possible, with a zero-bias feature retained even when $k_BT_K \ll \Delta$ ~\cite{Clerk2000,Cuevas2001,Tanaka2007,Domanski2008}~, most experiments carried out so far on N-QD-S and S-QD-S systems by various groups ~\cite{Buitelaar2002,Buizert2007,Yeyati2003,Sand-Jespersen2007,Graber2004,Deacon2010a,Deacon2010b,Kanai2010}~ do not observe Kondo features at this limit.

Here we suggest that the co-existence between the Kondo effect and superconductivity is a consequence of the 2-gap nature of NbSe$_2$\xspace.
Since we observe the Kondo zero-energy state even at fields far from critical fields of bulk 2H-NbSe$_2$\xspace, $H_{c2}^{\perp} \approx $ 4 T, $H_{c2}^{\parallel} \approx$ 17 T~\cite{Xi2016}~, the normal electrons must be associated with a superconducting band that turns normal even at such low fields.
As discussed previously, bulk NbSe$_2$\xspace has two superconducting gaps - one associated with niobium based bands around the $\Gamma$ and $K$ points and a second proximitized inner gap from the small selenium based bands around the $\Gamma$ point ~\cite{Dvir2018a,Kiss2007}~ whose signatures are observed at 1.25 mV and 0.6 mV respectively.

To elucidate the role of the Se-derived band in mediating the Kondo effect, we examine the correlation between the appearance of a zero-energy state and the second band in out-of-plane magnetic field (Figure 4). The out-of-plane magnetic field is chosen since the inner gap exhibits stronger response to its onset~\cite{Dvir2018b}~. Figure 4(a) shows the zero-bias conductance as a function of $B_\perp$ as it is swept in the positive direction (indicated by the red arrow). A complete dataset tracing the zero-bias conductance $G_0$ from negative to positive fields is shown in Supplementary Figure 1. The conductance reaches the lowest value of 0.16 $\mu$S at $B_\perp$ = -0.2 mT, where the tunneling spectrum reaches a sub-gap conductance typical of a hard-gap NbSe$_2$\xspace junction ~\cite{Dvir2018a}~. At this limit, there is no indication of a zero-bias conductance peak.

As $B_\perp$ is increased, $G_0$ increases in discrete steps of $\approx$ 1 $\mu$S, observed at 1.0 mT, 2.2 mT and 3.1 mT. As $B_\perp$ is further increased, the $G_0$ increase becomes continuous. Tunneling spectra are measured for fields that mark a discrete increase in $G_0$ (Figure 4(b)) and also at two higher fields. The second derivative of tunneling current as a function of magnetic field is plotted in Figure 4(c). The data in Figure 4(b)-(c) are offset along the $y$-axis for clarity. We observe that the first signature of a distinct zero-bias peak appears at a field of 3.1 mT and persists at higher fields. Simultaneously, the feature at 0.6 mV corresponding to the Se derived superconducting bands loses prominence until it is difficult to resolve from its background.

\begin{figure}[ht]
    \centering
    \includegraphics[width = 400 pt]{figure4_v11.png}
    \caption{(a) Differential conductance at zero-bias ($G_0$) as a function of $B_\perp$ at mT resolution. The red arrow indicates the sweep direction of the $B_\perp$; (b) sub-gap tunneling spectra at various field values indicated by the coloured dotted lines in (a); (c) second derivative of the tunneling current at field values indicated by dotted lines in (a); (d) zero-bias conductance peak height (left $y$-axis) and height of the feature at 0.6 mV in second derivative of tunneling current (right $y$-axis) as a function of $B_\perp$}
    \label{fig:fig4}
\end{figure}

Figure 4(d) shows the evolution of the zero-bias conductance peak height (left $y$-axis) and the height of the feature in the second derivative of tunneling current (right $y$-axis) corresponding to the inner Se superconducting band at 0.6 mV. An inverse correlation is evident from the data: The peak at 0.6 mV is the strongest at -0.2 mT where the conductance at zero-bias reaches the lowest value. With reduction in height of the feature in magnetic field, the zero-bias conductance feature grows in strength. The reduced height of the zero-bias conductance peak at higher fields can be attributed to the higher background conductance with magnetic field. Supplementary Figure 2 shows the individual sub-gap tunneling spectra shown in Figure 4(b) with their background conductance. The discrete jumps ($\approx 1\mu$S) in the zero-bias conductance with magnetic field can be associated with the discrete entry of vortices in the tunnel junction area of NbSe$_2$\xspace ~\cite{Dvir2018b}~. This correlation thus supports the interpretation that the Se band turning  normal could assist in the spin-flip co-tunneling event, as depicted in the schematic Figure 5(a).
Interestingly, we have shown in a previous work that in the weak coupling limit, a defect-bound QD favors coupling to the Se-band over coupling to the Nb-bands~\cite{Devidas2021}~.

\begin{figure}[ht]
    \centering
    \includegraphics[width = 400 pt]{figure5_v16.png}
    \caption{(a) Schematic of spin-flip co-tunneling event – (above) prohibited at zero applied magnetic field due to absence of states within the superconducting gap, (below) allowed via filled second band in NbSe$_2$\xspace in an applied magnetic field; (b) Sub-gap density of states of NbSe$_2$\xspace as a function of $B_{||}$, obtained via a large N expansion with a broadening $\delta$ = 0.04 meV. The inset shows the experimental differential tunneling conductance data for the same values of $B_{||}$.}
    \label{fig:fig5}
\end{figure}

To better understand the origin of the observed zero-bias conductance peak we
study a Kondo model in which we
treat the weakly coupled QD
as an impurity, with magnetic moment $\ssb$, placed at position $\bR$.
The Hamiltonian for the electrons in NbSe$_2$\xspace can
be written as $H=H_0+H_J$, where $H_0$ is the Hamiltonian
when no QD is present, and $H_J$ describes the coupling of the QD to the electrons in NbSe$_2$\xspace.
In the Kondo limit
$H_J= J \sum_{\sigma\sigma'} {c}^\dagger_{\bR\sigma}{\bm\tau}_{\sigma\sigma'}{c}_{\bR\sigma'}\cdot\ssb$, where
$J>0$ is the antiferromagnetic coupling between the QD and NbSe$_2$ electronic states, $\tau_i$ are the Pauli
matrices in spin space,
${c}^\dagger_{\bR\sigma}$ and ${c}_{\bR\sigma}$ are the creation and annihilation operators, respectively, for an electron at position $\bR$ with spin $\sigma$.
To treat the interaction between $\ssb$ and the electrons we use the large-N expansion~\cite{Read1983,Bickers1987}~.
The large-N expansion has been shown to give accurate results, in agreement with other methods such
as the numerical renormalization group~\cite{bulla2008a,gonzalez-buxton1998}~,
when the impurity is isolated and the Kondo problem has only one channel.
The electrons' density of states $\rho(\varepsilon)$ is the key property that determines the features of the Kondo effect.
Here, instead of assuming a constant $\rho(\varepsilon)$, or a prescribed form~\cite{Rossi2006,Mastrogiuseppe2014,Principi2015}~,
we use the density of states extracted for different values of the magnetic field
from spectroscopic measurements on few-layer superconducting NbSe$_2$\xspace~\cite{Dvir2018a}~.
One important aspect of using the experimentally obtained $\rho(\varepsilon)$ is that it allows us
to correctly capture the effect of $B$ on it, and therefore the unusual evolution
of the Kondo peak with $B$.
In the large-N expansion $\ssb$ is expressed in terms
of auxiliary creation (annihilation) fermionic operators
${f}^\dagger_\sigma$ (${f}_\sigma$)
satisfying the constraint
$n_f=\sum_\sigma {f}^\dagger_\sigma {f}_{\sigma}=1$,
with  $\sigma = 1,\ldots, N$, so that
$H_J = J \sum_{\kk,\kk',\sigma,\sigma'} {c}^\dagger_{\kk\sigma}{c}_{\kk'\sigma'}{f}^\dagger_{\sigma'}{f}_\sigma$, with
${c}^\dagger_{\kk\sigma}$ (${c}_{\kk\sigma}$) the creation (annihilation) operators for
an electron with momentum $\kk$ and spin $\sigma$.
In the remainder we assume $|\ssb|=1/2$ and therefore set $N=2$.
We decouple the quartic interaction term $H_J$ via the mean-field
$m\sim \sum_{\kk,\sigma} \langle  {\hat f}^\dagger_{\sigma} {\hat c}_{\kk\sigma} \rangle$.

The constrain $|\ssb|=1/2$ implies $n_f=1$ and this is enforced via a Lagrange multiplier, $\mu_f$,
which plays the role of the chemical potential of the $f$-electrons.
By minimizing the effective action corresponding to $H=H_0+H_J$, within the saddle-point approximation,
we obtain~\cite{Bickers1987}~
\begin{align}
 &\sum_{\sigma=\pm 1}\int_{-D}^{D}d\varepsilon n_F(\varepsilon)
 \frac{\rho(\varepsilon)(\varepsilon-\mu_f+h\sigma)}{(\varepsilon-\mu_f+h\sigma)^2+[\pi\rho(\varepsilon)m^2]^2}+\frac{2}{J}=0;
 \label{eq-m} \\
 &\sum_{\sigma=\pm 1}\int_{-D}^{D}d\varepsilon n_F(\varepsilon)
 \frac{\rho(\varepsilon)m^2}{(\varepsilon-\mu_f + h\sigma)^2 + [\pi\rho(\varepsilon)m^2]^2} -  1 = 0
 \label{eq-muf}
\end{align}
where $n_F(\varepsilon)$ is the Fermi function,
$D$ is the electrons' bandwidth, and
$h\approx (1/2)g_f\mu_B B$,  describes the Zeeman effect, with $g_f$ the effective g-factor for the QD
extracted from the experimentally observed Zeeman splitting.
Equations~\ceq{eq-m},~\ceq{eq-muf} can be inverted to find
$m^2$ and $\mu_f$ for a given value of $J$ and $B$.

From Eqs.~\ceq{eq-m},~\ceq{eq-muf} we obtain the values of $m$ and $\mu_f$ in the limit of $T\to 0$ (we set $T=6\times10^{-4}$~K).
For $\varepsilon\gg\Delta_1$, $\rho(\varepsilon)$, to very good approximation, can be taken to be constant.
Let $\rho_0$ be the value of $\rho(\varepsilon)$ for $\varepsilon\gg\Delta_1$.
The value of $\rho_0$ that enters Eqs.~\ceq{eq-m},~\ceq{eq-muf} depends on the effective spatial extension of the magnetic
impurity, i.e., in our case, of the defect-bound QD.
$\rho_0$ and $J$ are very difficult to estimate, both theoretically and experimentally.
On the other hand, $m$ and $\mu_f$, depend mostly on the product $\rho_0 J$,
and not on the separate values of $\rho_0$ and $J$. We studied the dependence of
the Kondo peak on the product of $\rho_0 J$ and found the general result that the peak grows with magnetic field for
$\rho_0 J\lesssim 0.2$ and decreases for $\rho_0 J\gtrsim 0.2$.
We find that the scaling of the Kondo peak with $B$ observed experimentally is best approximated
when $\rho_0 J\approx 0.14$. In the remainder we therefore set $\rho_0 J = 0.14$.

We obtain the electrons' Green's function at position $\bR$ and energy $\omega$, $G(\bR,\bR,\omega)$, renormalized by the coupling to the effective magnetic impurity, via the equations
\begin{align}
 & G(\bR,\bR,\omega) = G_0(\bR,\bR,\omega) + m^2 G_0(\bR,\bR,\omega) F(\omega) G_0(\bR,\bR,\omega),
 \label{eq-G} \\
 & F(\omega)=[(\omega+i\delta-\mu_f-m^2 G_0(\bR,\bR,\omega))\tau_0-(1/2)g\mu_B B\tau_z]^{-1},
 \label{eq-F}
\end{align}
where
$G_0(\bR,\bR,\omega)=-i(\pi/2)\rho(\omega)\tau_0$,
is the electrons' bare local Green's function,
$F(\omega)$
is the Green's function for the effective fermionic degrees of freedom $f$, and $\delta$ is
a small broadening that we introduce to take into account the effect on the spectra of thermal fluctuations and non-magnetic disorder.
The local DOS at the position of the QD, at energy $\omega$, $\rho_{\rm tot}(\omega)$ is then obtained as
$\rho_{\rm tot}(\omega)=-{\Im}[{\rm Tr}(G)]/\pi$.

Figure~\ref{fig:fig5}~(b) shows  $\rho_{\rm tot}(\omega)$ for different values of $B$ obtained via the large-N expansion assuming
$\delta=0.04 ~\text{meV}$. We see that for $B=0$, given the vanishing of the DOS as $\varepsilon\to 0$ for superconducting NbSe$_2$\xspace, no peak for $\varepsilon\approx 0$ is present in $\rho_{\rm tot}(\omega)$, indicating the absence of Kondo screening.
NRG results~\cite{Ingersent1996}~ show that when $N(\varepsilon)\sim |\varepsilon|^a$ with a $>$ 1/2, and perfect particle-hole symmetry is present, no Kondo effect can be realized.
In superconducting NbSe$_2$\xspace the presence of a magnetic field $B$ induces a finite density of states at low energies.
This makes possible the establishment of a Kondo cloud, and, when $B$ is not too large,
dominates over the suppression of $T_K$ due to the polarization of $\ssb$ induced by the magnetic field.
For $B\gtrsim 0.25$~T a peak develops: the softening, at low energies, of the superconducting
gap allows the establishment of Kondo screening. At larger $B$ the Kondo peak splits, as expected~\cite{Hewson1993}~.

\section{Summary}
Our work shows that defect-bound QDs can be used as platforms to probe the interplay between the energy scales associated with the Kondo effect and with superconductivity. This opens a number of interesting future possibilities. First, as discussed above, NbSe$_2$\xspace is a 2-gap superconductor - offering a richer phase space involving both order parameters. Interestingly, NbSe$_2$\xspace is expected to develop a triplet order parameter at high $B_\parallel$~\cite{mockli2020}~. This could give rise to a field-driven transition in the type of Kondo effect. Finally, it is interesting to consider coupling QDs to unconventional SCs such as FeTe$_{0.55}$Se$_{0.45}$~\cite{Zalic2019}~, which can be produced using exfoliation, or even
to the superconducting states appearing in twisted graphene systems~\cite{Cao2018a}~.

\section{Methods}
The tunneling measurements are performed using standard low frequency lockin techniques in a BlueFors dilution cryostat with a base temperature of 20 mK. An AC excitation of 30$\mu V$ is applied across the device between the tunneling electrodes and the ohmic drain electrode using a Zurich instruments MFLI digital lock-in amplifier (LIA). The superconductor is tuned in and out of the superconducting gap by applying a DC bias voltage ($V_{SD}$), also obtained from the LIA. The tunneling current reaching the drain is fed into a FEMTO current amplifier. The output of the current amplifier is then fed into a Keithley 2000 digital multimeter to measure $I_{DC}$ and into the LIA input to measure the differential conductance $dI/dV$ as a function of $V_{SD}$. The magnetic field measurements are performed with the 2 axis 9T-3T (Z-Y) superconducting magnet attached to the dilution. The magnet is powered by the commercial power supply from American Magnetic Inc. (AMI), with a resolution of 1 mT (supplied). For resolutions smaller than that, the magnet is powered by a Keithley 2400 SMU in the current source mode. The applied field is evaluated from the current to field ratio provided by AMI.

\section{Acknowledgements}
Devices for this study were fabricated at the Center for Nanoscience and Nanotechnology, The Hebrew University.
Funding for this work was provided by a European Research Council Starting Grant (No. 637298, TUNNEL), Israel Science Foundation grant 861/19, and BSF grant 2016320.
E.R.  thanks the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611, where part of this work was conducted, and Kevin Ingersent for very helpful discussions.


\section{Introduction}

Athletes benefit from feedback from professional evaluations to maximize performance and minimize injury risk.
Often, these evaluations require athletes to perform key movements that facilitate evaluation, e.g. drop or squat jumps, with the use of specialized equipment like force plates. However, not all athletes have access to this equipment, nor to trained professionals to help them understand the variables measured by these movements and how they should adjust their training in response. Here, we establish that key variables from these jumps can be estimated using from videos of athletes performing both squat and drop jumps.

We use force plate data as the ground truth for our variables. The accuracy of force plates has been well established. Garcia-Ramos \etal \cite{garcia2017reliability} show that the force plate variables like velocity and power are highly reliable when evaluating a countermovement jump for an athlete's performance. Velocity is an important variable when determining posture \cite{duarte2010revision} or to represent the maximum lower-limb capacity of an athlete \cite{samozino2014force}, which is helpful in creating individual and prioritized training profiles for athletes.

\begin{figure}[htp!]
\begin{center}
   \includegraphics[width=0.45\textwidth]{arch.drawio.pdf}
\end{center}
   \caption{Overview of our pipeline. The Swin Transformer block takes as input the video frames and gives out flattened features further processed through two dense layers for an estimated velocity score.}
\label{fig:arch}
\end{figure}

The usefulness of force plates is limited to laboratory settings \cite{johnson2018predicting}, but it is substantially more feasible to record a video outside of such settings and utilize a computer vision model to estimate these key variables. However, for this to be realistic, the computer vision model has to be sufficiently reliable. Prior work has examined the feasibility of identifying when an athlete incorrectly performs a drop or squat jump \cite{roygaga2022ape}, but thus far no works have estimated the key variables themselves.

Prior work has also focused on alternate solutions for athlete performance evaluation, useful in the absence of standard force plate data. Johnson \etal~\cite{johnson2018predicting} provide a method to analyze athlete ground reaction forces and moments (GRF/Ms) in real time using eight marker motion capture trajectories, allowing for a field-based study in identifying early injury risk; they \cite{johnson2019field} also use the motion capture data to estimate Knee Joint Movements (KJM). Goldacre \etal \cite{goldacre2021predicting} utilized preprocessed marker motion capture data and compared a Least Squares Estimator method for estimating GRF with a convolutional neural network (CNN) estimator and found the CNN-best method was more reliable and consistent. All of these works showcase the feasibility of estimating a key metric, GRF, but none have considered other essential metrics like velocity. Further, none of these works have attempted to estimate variables using raw video data. Here, we train a custom Swin Transformer architecture to estimate key metrics like velocity using the dataset from Blanchard \etal \cite{blanchard2019keep}.

We propose an automatic force plate variable estimation tool that uses athlete motion videos --- performing evaluative jumps such as countermovement and drop jumps ---
as input, which can be run on any handy smart device such as a mobile phone. The tool is trained on the videos directly with minimal preprocessing steps, allowing faster inference during its execution. In the absence of a force plate, our method can provide a good approximation of the desired force plate variables, which can further be used to investigate the performance of an athlete.

In summary, this work has several novel contributions:
\begin{enumerate}
    \item This is the first work to showcase that key variables for assessing athlete performance can be estimated using video alone.
    \item We create a custom variant on the Swin Transformer architecture (see Figure \ref{fig:arch}) and achieve state-of-the-art performance across 86 athletes.
    \item Our solution is generalized across two different kinds of jumps: the squat jump and the drop jump. We did not create separate models for each kind of jump.
\end{enumerate}

\section{Methods}

\subsection{Data}

Blanchard \etal \cite{blanchard2019keep} collected a dataset which consisted of 582 RGB videos of 89 athletes performing multiple countermovement and drop jumps~\cite{blanchard2019keep}. Each jump was recorded from three different angles: a left, right, and center perspective. Jumps were performed on force plate data, which recorded a multitude of variables associated with the jump (181 for countermovement jumps and 58 for drop jumps). Only 86 athletes' force plate data (545 videos) were available, so we only utilized this subset of the data.

We divided the dataset into three separate folds. No participant from one fold was present in any other fold. The first two folds included 29 unique participants each, and the third fold included the remaining 28 participants.
There were 180, 178, and 161 videos in each fold, respectively.

We focused on training a model to estimate concentric peak velocity. Figure~\ref{fig:Fold Graph} displays the distribution of concentric peak velocity (hereafter, velocity) across the three folds of our dataset.
In our experiments, models used for evaluation were trained on each individual view, with an additional model trained on the combined input of the three views.
\begin{figure}[htp!]
\begin{center}
   \includegraphics[width=1\linewidth,trim=4 4 4 4,clip] {COMBFIG.pdf}
\end{center}
   \caption{Distribution of velocity across data folds.}
\label{fig:Fold Graph}
\end{figure}

\subsection{Model Architectures}

Estimating velocity, as well as other metrics from the force plate data, is a regression problem. We experimented with three different methods to evaluate the feasibility of estimating velocity.

\textbf{Baseline Models} We evaluated two different baseline methods for estimating the participant's velocity using video data in order to measure athletic performance without a force plate setup. First, we used a ResNet-50 architecture to extract features from all video frames. For our first baseline, we trained a support vector machine (SVM) using these features. For our second, we added a regression dense layer to the ResNet-50.

In the first baseline method, we use the SVM regression method for velocity estimation, with the kernel set to 'rbf', degree to 3, coef0 to 0.0, tol to 0.0001, C to 1.0 and  epsilon to 0.01. The second baseline method uses a dense layer for final score estimation.

\textbf{Swin Transformer Architecture} We hypothesized a transformer architecture~\cite{dosovitskiy2020image} would be better suited for video data than our baseline methods. We modified the \textit{Swin Transformer}~\cite{liu2021swin} to extract features from video frames and utilized a dense layer to predict velocity. The first layer of the Swin Transformer divided the image into smaller patches with equal window sizes using the window multi-attention head (W-MSA), as shown in Equation~\ref{eq1}. The second layer repeated this, but shifted the starting position using shifted window multi-attention head (SW-MSA), as shown in Equation~\ref{eq2}. This change in how window patches are made across layers produced cross-connections that let the model identify the relevant image features. As illustrated in Equations~\ref{eq3} and~\ref{eq4}, self attention was conducted in the local window patches before the window patches were sent via a multi-layer perceptron (MLP). Finally, the adjacent patches were combined. The two successive layers were part of a Swin Transformer block; we used four such Swin Transformer blocks to generate the features with respect to each video frame.

\begin{equation}
\hat{z}^{l}= \mathrm{W}\\-\mathrm{MSA} \left( \mathrm{LN}\left( \mathrm{z}^{l-1}\right)\right) + \mathrm{z}^{l-1}
\label{eq1}
\end{equation}

\begin{equation}
\hat{z}^{l+1}= \mathrm{SW-MSA} \left( \mathrm{LN}\left( \mathrm{z}^{l}\right)\right) + \mathrm{z}^{l}
\label{eq2}
\end{equation}

\begin{equation}
z^{l}= \mathrm{MLP} \left( \mathrm{LN}\left( \mathrm{\hat{z}}^{l}\right)\right) + \mathrm{\hat{z}}^{l}
\label{eq3}
\end{equation}

\begin{equation}
z^{l+1}= \mathrm{MLP} \left( \mathrm{LN}\left( \mathrm{\hat{z}}^{l+1}\right)\right) + \mathrm{\hat{z}}^{l+1}
\label{eq4}
\end{equation}

where $\hat{z}^{l}$ are the output features of window-based multihead self-attention layers ($\mathrm{W}-\mathrm{MSA}$) and shifted window-based multihead self-attention layers ($\mathrm{SW}-\mathrm{MSA}$) for block $l$ respectively.

The architecture with the Swin Transformer backbone is as shown in Figure~\ref{fig:arch}. The inputs to the model are video frames from each view with respect to a single participant. Each of the video frames are processed through the Swin Transformer Block~\cite{liu2021swin}.

The mean of the frames from each view after flattening these features are calculated. Each view is passed through a dense layer before concatenating them all together. The combined features are then processed through another dense layer to obtain the final features.

\subsection{Training and Evaluation Details}

Velocity scores followed a Gaussian distribution, with the majority of data ranging between 0.4 and 0.6, as seen in Figure \ref{fig:Fold Graph}. We exploited this distribution during training by using a balanced sampler to create mini batches, thus allowing the model to accurately learn the embedding space of jumps across all velocities.

All of our models were trained on a Titan V GPU with an Intel(R) Xeon(R) E5-1650 v4 processor and 16GB RAM. We used L1 loss for training and kept a fixed learning rate of 1e-3. The models were trained for 100 epochs using the Adam Optimizer.

All models were trained and evaluated following a 3-fold cross validation paradigm. No participants present in the training data were present in the test data.

We evaluated models using two metrics: Pearson's R and Mean Absolute Error (MAE). All reported results are the mean and standard deviation of each metric across all folds.

\section{Results}

\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Our baseline when compared to SVM and Dense Layers for estimation on the Video and Force plate data.}
\label{table:model_baseline_compare}
\begin{sc}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| >{\raggedright\arraybackslash}p{2.7cm} | >{\raggedright\arraybackslash}p{2cm} | >{\raggedright\arraybackslash}p{2cm} |}
 \hline
 {\bfseries Method} & {\bfseries MAE ($\downarrow$)} & {\bfseries R  ($\uparrow$)}\\
 \hline
 {Baseline method with SVM} & {0.19 \textpm 0.017} & {0.12 \textpm 0.06}\\
 \hline
 {Baseline method with Dense Layer} & {0.28 \textpm 0.017 } & {0.16 \textpm 0.05 }\\
 \hline
 \hline
 \textbf{Our method (Combined view)} & {\bfseries 0.21 \textpm 0.015} & {\bfseries 0.71 \textpm 0.06}\\
 \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vspace{-0.2cm}
\end{table}

Table~\ref{table:model_baseline_compare} shows the performance of our baseline estimators and our transformer method. Although these methods appear similar when considering only MAE, our transformer-based model performs substantially better when evaluated by Pearson's R. In Figure 3, we showcase how our estimated velocity compares with the actual velocity for each fold --- our models slightly underestimated velocity, but they were extremely consistent with their estimations across folds. The low R values of the baseline methods indicate the estimations from those models were extremely inconsistent.

\begin{figure}[htp!]
     \centering
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold1.pdf}
         \caption{Fold 1}
         \label{fig:V1}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold2.pdf}
         \caption{Fold 2}
         \label{fig:V2}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold3.pdf}
         \caption{Fold 3}
         \label{fig:V3}
     \end{subfigure}
        \caption{Actual vs estimated velocity scores for the Transformer pipeline trained on a combined view input.}
        \label{fig:Fold Graph1}
\end{figure}

\subsection{Ablation}

We investigated how influential each component of the data was for our method's prediction. Table~\ref{table:model_view_compare} showcases the results of this experiment --- video from the center perspective was more important for accurately estimating velocity, exhibiting an R value 0.21 higher than either the left or right perspective (the left and right perspectives exhibited similar performances). This contrasts the findings of \cite{roygaga2022ape}, which found the left and right perspectives were more influential for estimating when a jump was performed incorrectly. This indicates that different perspectives may be essential to fully estimating all components of a jump. Further, data from all three perspectives culminated in the best model --- combining all perspectives increased the R value by 0.15 over the center view alone.

\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Experiments to determine the importance of combining videos captured from different views. The transformer pipeline trained on the combination of the three views improves the correlation score (R), which shows that the model is able to distinguish between a low and a high velocity jump when provided with video input from multiple angles. But the model trained on the Left view provides the least error in estimation, followed by the Right view, which shows that the model can still learn accurate estimations from videos from a single view but with a much lower capacity to distinguish between the higher and lower velocity jumps.}
\label{table:model_view_compare}
\begin{sc}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| >{\raggedright\arraybackslash}p{2.7cm} | >{\raggedright\arraybackslash}p{2cm} | >{\raggedright\arraybackslash}p{2cm} |}
 \hline
 {\bfseries Video capture view} & {\bfseries MAE ($\downarrow$)} & {\bfseries R  ($\uparrow$)}\\
 \hline
 {Left} & {\bfseries 0.16 \textpm 0.016} & {0.35 \textpm 0.06}\\
 \hline
 {Right} & {0.20 \textpm 0.014} & {0.34 \textpm 0.04}\\
 \hline
 {Center} & {0.19 \textpm 0.014} & {0.56 \textpm 0.06}\\
 \hline
 {\textbf{Combined View}} & {0.21 \textpm 0.015} & {\bfseries 0.71 \textpm 0.06}\\
 \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vspace{-0.2cm}
\end{table}

\section{Conclusion and Future Work}

Velocity is a key evaluation metric for evaluating an athlete's performance. To our knowledge, this is the first work that estimates velocity from RGB video alone. Our novel solution is generalized across participants and two kinds of evaluative jumps: the squat jump and the drop jump. This work also provides insights into what kinds of data are important for estimating key metrics like velocity. In the future, we plan to expand our architecture to estimate a multitude of other essential variables for athlete performance following a multi-task regression paradigm.

{\small
\bibliographystyle{ieee_fullname}

\section{Introduction}

Athletes benefit from feedback from professional evaluations to maximize performance and minimize injury risk.
Often, these evaluations require athletes to perform key movements that facilitate evaluation, e.g. drop or squat jumps, with the use of specialized equipment like force plates. However, not all athletes have access to this equipment, nor to trained professionals to help them understand the variables measured by these movements and how they should adjust their training in response. Here, we establish that key variables from these jumps can be estimated using from videos of athletes performing both squat and drop jumps.

We use force plate data as the ground truth for our variables. The accuracy of force plates has been well established. Garcia-Ramos \etal \cite{garcia2017reliability} show that the force plate variables like velocity and power are highly reliable when evaluating a countermovement jump for an athlete's performance. Velocity is an important variable when determining posture \cite{duarte2010revision} or to represent the maximum lower-limb capacity of an athlete \cite{samozino2014force}, which is helpful in creating individual and prioritized training profiles for athletes.

\begin{figure}[htp!]
\begin{center}
   \includegraphics[width=0.45\textwidth]{arch.drawio.pdf}
\end{center}
   \caption{Overview of our pipeline. The Swin Transformer block takes as input the video frames and gives out flattened features further processed through two dense layers for an estimated velocity score.}
\label{fig:arch}
\end{figure}

The usefulness of force plates is limited to laboratory settings \cite{johnson2018predicting}, but it is substantially more feasible to record a video outside of such settings and utilize a computer vision model to estimate these key variables. However, for this to be realistic, the computer vision model has to be sufficiently reliable. Prior work has examined the feasibility of identifying when an athlete incorrectly performs a drop or squat jump \cite{roygaga2022ape}, but thus far no works have estimated the key variables themselves.

Prior work has also focused on alternate solutions for athlete performance evaluation, useful in the absence of standard force plate data. Johnson \etal~\cite{johnson2018predicting} provide a method to analyze athlete ground reaction forces and moments (GRF/Ms) in real time using eight marker motion capture trajectories, allowing for a field-based study in identifying early injury risk; they \cite{johnson2019field} also use the motion capture data to estimate Knee Joint Movements (KJM). Goldacre \etal \cite{goldacre2021predicting} utilized preprocessed marker motion capture data and compared a Least Squares Estimator method for estimating GRF with a convolutional neural network (CNN) estimator and found the CNN-best method was more reliable and consistent. All of these works showcase the feasibility of estimating a key metric, GRF, but none have considered other essential metrics like velocity. Further, none of these works have attempted to estimate variables using raw video data. Here, we train a custom Swin Transformer architecture to estimate key metrics like velocity using the dataset from Blanchard \etal \cite{blanchard2019keep}.

We propose an automatic force plate variable estimation tool that uses athlete motion videos --- performing evaluative jumps such as countermovement and drop jumps ---
as input, which can be run on any handy smart device such as a mobile phone. The tool is trained on the videos directly with minimal preprocessing steps, allowing faster inference during its execution. In the absence of a force plate, our method can provide a good approximation of the desired force plate variables, which can further be used to investigate the performance of an athlete.

In summary, this work has several novel contributions:
\begin{enumerate}
    \item This is the first work to showcase that key variables for assessing athlete performance can be estimated using video alone.
    \item We create a custom variant on the Swin Transformer architecture (see Figure \ref{fig:arch}) and achieve state-of-the-art performance across 86 athletes.
    \item Our solution is generalized across two different kinds of jumps: the squat jump and the drop jump. We did not create separate models for each kind of jump.
\end{enumerate}

\section{Methods}

\subsection{Data}

Blanchard \etal \cite{blanchard2019keep} collected a dataset which consisted of 582 RGB videos of 89 athletes performing multiple countermovement and drop jumps~\cite{blanchard2019keep}. Each jump was recorded from three different angles: a left, right, and center perspective. Jumps were performed on force plate data, which recorded a multitude of variables associated with the jump (181 for countermovement jumps and 58 for drop jumps). Only 86 athletes' force plate data (545 videos) were available, so we only utilized this subset of the data.

We divided the dataset into three separate folds. No participant from one fold was present in any other fold. The first two folds included 29 unique participants each, and the third fold included the remaining 28 participants.
There were 180, 178, and 161 videos in each fold, respectively.

We focused on training a model to estimate concentric peak velocity. Figure~\ref{fig:Fold Graph} displays the distribution of concentric peak velocity (hereafter, velocity) across the three folds of our dataset.
In our experiments, models used for evaluation were trained on each individual view, with an additional model trained on the combined input of the three views.
\begin{figure}[htp!]
\begin{center}
   \includegraphics[width=1\linewidth,trim=4 4 4 4,clip] {COMBFIG.pdf}
\end{center}
   \caption{Distribution of velocity across data folds.}
\label{fig:Fold Graph}
\end{figure}

\subsection{Model Architectures}

Estimating velocity, as well as other metrics from the force plate data, is a regression problem. We experimented with three different methods to evaluate the feasibility of estimating velocity.

\textbf{Baseline Models} We evaluated two different baseline methods for estimating the participant's velocity using video data in order to measure athletic performance without a force plate setup. First, we used a ResNet-50 architecture to extract features from all video frames. For our first baseline, we trained a support vector machine (SVM) using these features. For our second, we added a regression dense layer to the ResNet-50.

In the first baseline method, we use the SVM regression method for velocity estimation, with the kernel set to 'rbf', degree to 3, coef0 to 0.0, tol to 0.0001, C to 1.0 and  epsilon to 0.01. The second baseline method uses a dense layer for final score estimation.

\textbf{Swin Transformer Architecture} We hypothesized a transformer architecture~\cite{dosovitskiy2020image} would be better suited for video data than our baseline methods. We modified the \textit{Swin Transformer}~\cite{liu2021swin} to extract features from video frames and utilized a dense layer to predict velocity. The first layer of the Swin Transformer divided the image into smaller patches with equal window sizes using the window multi-attention head (W-MSA), as shown in Equation~\ref{eq1}. The second layer repeated this, but shifted the starting position using shifted window multi-attention head (SW-MSA), as shown in Equation~\ref{eq2}. This change in how window patches are made across layers produced cross-connections that let the model identify the relevant image features. As illustrated in Equations~\ref{eq3} and~\ref{eq4}, self attention was conducted in the local window patches before the window patches were sent via a multi-layer perceptron (MLP). Finally, the adjacent patches were combined. The two successive layers were part of a Swin Transformer block; we used four such Swin Transformer blocks to generate the features with respect to each video frame.

\begin{equation}
\hat{z}^{l}= \mathrm{W}\\-\mathrm{MSA} \left( \mathrm{LN}\left( \mathrm{z}^{l-1}\right)\right) + \mathrm{z}^{l-1}
\label{eq1}
\end{equation}

\begin{equation}
\hat{z}^{l+1}= \mathrm{SW-MSA} \left( \mathrm{LN}\left( \mathrm{z}^{l}\right)\right) + \mathrm{z}^{l}
\label{eq2}
\end{equation}

\begin{equation}
z^{l}= \mathrm{MLP} \left( \mathrm{LN}\left( \mathrm{\hat{z}}^{l}\right)\right) + \mathrm{\hat{z}}^{l}
\label{eq3}
\end{equation}

\begin{equation}
z^{l+1}= \mathrm{MLP} \left( \mathrm{LN}\left( \mathrm{\hat{z}}^{l+1}\right)\right) + \mathrm{\hat{z}}^{l+1}
\label{eq4}
\end{equation}

where $\hat{z}^{l}$ are the output features of window-based multihead self-attention layers ($\mathrm{W}-\mathrm{MSA}$) and shifted window-based multihead self-attention layers ($\mathrm{SW}-\mathrm{MSA}$) for block $l$ respectively.

The architecture with the Swin Transformer backbone is as shown in Figure~\ref{fig:arch}. The inputs to the model are video frames from each view with respect to a single participant. Each of the video frames are processed through the Swin Transformer Block~\cite{liu2021swin}.

The mean of the frames from each view after flattening these features are calculated. Each view is passed through a dense layer before concatenating them all together. The combined features are then processed through another dense layer to obtain the final features.

\subsection{Training and Evaluation Details}

Velocity scores followed a Gaussian distribution, with the majority of data ranging between 0.4 and 0.6, as seen in Figure \ref{fig:Fold Graph}. We exploited this distribution during training by using a balanced sampler to create mini batches, thus allowing the model to accurately learn the embedding space of jumps across all velocities.

All of our models were trained on a Titan V GPU with an Intel(R) Xeon(R) E5-1650 v4 processor and 16GB RAM. We used L1 loss for training and kept a fixed learning rate of 1e-3. The models were trained for 100 epochs using the Adam Optimizer.

All models were trained and evaluated following a 3-fold cross validation paradigm. No participants present in the training data were present in the test data.

We evaluated models using two metrics: Pearson's R and Mean Absolute Error (MAE). All reported results are the mean and standard deviation of each metric across all folds.

\section{Results}

\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Our baseline when compared to SVM and Dense Layers for estimation on the Video and Force plate data.}
\label{table:model_baseline_compare}
\begin{sc}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| >{\raggedright\arraybackslash}p{2.7cm} | >{\raggedright\arraybackslash}p{2cm} | >{\raggedright\arraybackslash}p{2cm} |}
 \hline
 {\bfseries Method} & {\bfseries MAE ($\downarrow$)} & {\bfseries R  ($\uparrow$)}\\
 \hline
 {Baseline method with SVM} & {0.19 \textpm 0.017} & {0.12 \textpm 0.06}\\
 \hline
 {Baseline method with Dense Layer} & {0.28 \textpm 0.017 } & {0.16 \textpm 0.05 }\\
 \hline
 \hline
 \textbf{Our method (Combined view)} & {\bfseries 0.21 \textpm 0.015} & {\bfseries 0.71 \textpm 0.06}\\
 \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vspace{-0.2cm}
\end{table}

Table~\ref{table:model_baseline_compare} shows the performance of our baseline estimators and our transformer method. Although these methods appear similar when considering only MAE, our transformer-based model performs substantially better when evaluated by Pearson's R. In Figure 3, we showcase how our estimated velocity compares with the actual velocity for each fold --- our models slightly underestimated velocity, but they were extremely consistent with their estimations across folds. The low R values of the baseline methods indicate the estimations from those models were extremely inconsistent.

\begin{figure}[htp!]
     \centering
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold1.pdf}
         \caption{Fold 1}
         \label{fig:V1}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold2.pdf}
         \caption{Fold 2}
         \label{fig:V2}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.3\textwidth}
         \centering
         \includegraphics[width=\textwidth]{Vel_fold3.pdf}
         \caption{Fold 3}
         \label{fig:V3}
     \end{subfigure}
        \caption{Actual vs estimated velocity scores for the Transformer pipeline trained on a combined view input.}
        \label{fig:Fold Graph1}
\end{figure}

\subsection{Ablation}

We investigated how influential each component of the data was for our method's prediction. Table~\ref{table:model_view_compare} showcases the results of this experiment --- video from the center perspective was more important for accurately estimating velocity, exhibiting an R value 0.21 higher than either the left or right perspective (the left and right perspectives exhibited similar performances). This contrasts the findings of \cite{roygaga2022ape}, which found the left and right perspectives were more influential for estimating when a jump was performed incorrectly. This indicates that different perspectives may be essential to fully estimating all components of a jump. Further, data from all three perspectives culminated in the best model --- combining all perspectives increased the R value by 0.15 over the center view alone.

\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Experiments to determine the importance of combining videos captured from different views. The transformer pipeline trained on the combination of the three views improves the correlation score (R), which shows that the model is able to distinguish between a low and a high velocity jump when provided with video input from multiple angles. But the model trained on the Left view provides the least error in estimation, followed by the Right view, which shows that the model can still learn accurate estimations from videos from a single view but with a much lower capacity to distinguish between the higher and lower velocity jumps.}
\label{table:model_view_compare}
\begin{sc}
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| >{\raggedright\arraybackslash}p{2.7cm} | >{\raggedright\arraybackslash}p{2cm} | >{\raggedright\arraybackslash}p{2cm} |}
 \hline
 {\bfseries Video capture view} & {\bfseries MAE ($\downarrow$)} & {\bfseries R  ($\uparrow$)}\\
 \hline
 {Left} & {\bfseries 0.16 \textpm 0.016} & {0.35 \textpm 0.06}\\
 \hline
 {Right} & {0.20 \textpm 0.014} & {0.34 \textpm 0.04}\\
 \hline
 {Center} & {0.19 \textpm 0.014} & {0.56 \textpm 0.06}\\
 \hline
 {\textbf{Combined View}} & {0.21 \textpm 0.015} & {\bfseries 0.71 \textpm 0.06}\\
 \hline
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vspace{-0.2cm}
\end{table}

\section{Conclusion and Future Work}

Velocity is a key evaluation metric for evaluating an athlete's performance. To our knowledge, this is the first work that estimates velocity from RGB video alone. Our novel solution is generalized across participants and two kinds of evaluative jumps: the squat jump and the drop jump. This work also provides insights into what kinds of data are important for estimating key metrics like velocity. In the future, we plan to expand our architecture to estimate a multitude of other essential variables for athlete performance following a multi-task regression paradigm.

{\small
\bibliographystyle{ieee_fullname}