\section{Introduction}
    The major breakthrough of the last 40 years in the field of superconductivity has been the discovery of several families of unconventional superconductors: heavy-fermions, organics, high-Tc cuprates, iron pnictides in chronological order. All are controlled by pairing mechanisms dominated by purely electronic interactions like magnetism, instead of the electron-phonon interaction.
    These new pairing mechanism also lead to new superconducting states, with different possible spin-states (spin-singlet/even parity or spin-triplet/odd parity), and corresponding orbital states, dubbed \dwave, \pwave, \fwave.
    Spin-triplet superconductivity is presently highly sought-after in the context of quantum engineering, as the required starting point to build robust topologically protected qubits.
    Artificial heterostructures have been proposed to create this rare state from conventional s-wave superconductors \cite{FuPRL2008}.
    However, experimental success following this line remains controversial.
    Hence bulk spin-triplet superconductors with strong spin-orbit coupling is an appealing alternative.

    Besides superfluid $^3$He \cite{Leggett1975}, a perfect analogue of a neutral spin-triplet \pwave superconductor below 1~mK, most candidates for bulk spin-triplet superconductivity are found among uranium-based heavy-fermion systems. The most famous is undoubtedly UPt$_3$, which is one of the rare systems  showing not only one superconducting phase, but three superconducting phases differing by their symmetries \cite{JoyntRMP2002,Sauls1994}.
    Another major series of systems is that of UGe$_2$, URhGe and UCoGe, where the bulk coexistence of ferromagnetism and superconductivity leaves little doubt that they are also  spin-triplet superconductors.
    They all display in addition uncommon superconducting phase diagrams, with a re-entrance or reinforcement of superconductivity under large magnetic fields  \cite{AokiRevJPSJ2019}.

    Naturally, these uranium heavy-fermion superconductors are also central in modern condensed matter physics, as potential topological superconductors.
    For example, UPt$_3$ is also considered as a prime candidate for topological chiral spin-triplet superconductivity, like URu$_2$Si$_2$ is a prime candidate for topological chiral spin-singlet (d-wave) superconductivity \cite{GhoshJPhysCondMat2020}.

    Nevertheless, spin-triplet superconductivity remains a rare state of matter, with no clear understanding of the main reasons driving its appearance.
    In this context, the discovery of superconductivity in \UTe above 1.5~K back in 2018 \cite{RanScience2019} drove much enthusiasm in the condensed matter physics community.
    \UTe is a stunning system for almost all of its properties, starting with the fact that it is metallic only thanks to strong electronic correlations \cite{Xu2019, IshizukaPRL2019, Shick2021}. Its superconductivity shows extraordinary robustness to very large magnetic fields, \cite{RanScience2019, Aoki2019, Knebel2019, RanNatPhys2019}, a feature traditionally associated to the ferromagnetic superconductors, and requiring at least spin-triplet superconductivity.
    Still \UTe is a paramagnetic heavy-fermion and we know now, with no evidence for being close to a ferromagnetic instability \cite{AokiJPCM2022}.
    Among the hardest yet most important questions is to understand what could drive such a system toward a spin-triplet superconducting ground-state ?
    Having better hints of where to look for or how to build spin-triplet superconductors would be of interest far beyond the community of quantum materials.
    This work brings unprecedented elements directly related to this issue, supporting notably that magnetic field  is likely driving \UTe from a spin-triplet to a spin-singlet ground state.

    Understanding how requires a few more explanations on the specificity of the system.
    Indeed, owing to these exceptional properties: field-reinforcement of superconductivity for field applied along the hard magnetisation \baxis \cite{RanScience2019, Knebel2019}, field-induced superconducting phase (between 40 and 60~T) for field applied in particular directions between the $b$ and \caxis \cite{RanNatPhys2019, Knafo2020},
    \UTe was immediately acknowledged as a spin-triplet superconductor, probably topological.

    Up to now, evidence for spin-triplet pairing in \UTe, comes from NMR Knight-shift measurements, recently performed along all crystallographic axes \cite{FujibayashiJPSJ2022}, and from the strong violation of the paramagnetic limit by the superconducting upper critical field \Hc \cite{RanScience2019, Knebel2019, RanNatPhys2019, Knafo2021, RanNPJQM2021}.
    Analysis of \Hc is however not straightforward in \UTe, because the field-reinforced phase observed for magnetic fields applied along the hard \baxis strongly suggests a field-enhanced pairing strength \cite{RanScience2019, Knebel2019}.
    Such a field-dependent pairing is difficult to model theoretically, and makes the analysis of \Hc more complex, because it opens other routes for a violation of the paramagnetic limit than spin-triplet pairing \cite{AokiJPCM2022}.
    Moreover, even the simplest question of whether or not the field-reinforced phase has a different symmetry than the low-field phase was up to now not settled.
    From the very beginning, the existence of this field-reinforced superconductivity inspired theoretical models for different superconducting order parameters between the low and high field phases.
    Yet experimentally no phase transition had been observed between these two "phases".

    Nonetheless, high pressure experiments demonstrate unambiguously that \UTe  does belong to the very select class of unconventional superconductors displaying transitions between different superconducting phases
    \cite{Braithwaite2019,Lin2020,Aoki2020,ThomasPRB2021}.
    Some samples even show a double transition at ambient pressure in zero field, which could explain how this low symmetry system (\UTe is orthorhombic) could also be chiral.
    Indeed, the observation of time reversal symmetry breaking below \Tc by polar Kerr-effect measurements \cite{Hayes2021}, or the asymmetric spectrum observed with STM spectroscopy \cite{Jiao2020} suggest such a chiral state.
    Actually, there is now more and more evidence that the double transition at ambient pressure in zero field is not an intrinsic effect \cite{AokiJPCM2022, ThomasPRB2021, Rosa2022}.

    This rich superconducting state inspired several theoretical work, proposing different scenarios for the possible symmetry states under magnetic field and pressure \cite{IshizukaPRL2019, Shishidou2021, Ishizuka2021}.
    Moreover, these works also address the already announced deeper question, of why this system is spin-triplet ?
    This question is particularly acute in \UTe, because contrary to initial expectations, no ferromagnetic fluctuations have yet been detected.
     Ferromagnetic spin-fluctuations, natural in the ferromagnetic superconductors \cite{AokiRevJPSJ2019} and necessary to explain the stability of the $A$-phase of superfluid $^3$He \cite{Leggett1975}, are natural candidates for the pairing mechanism of spin-triplet superconductors.
    By contrast in \UTe, microscopic neutron measurements have only revealed incommensurate magnetic fluctuations \cite{Duan2020, Knafo2021, ButchNPJ2022}, and a resonance at finite $\bm{Q}$ below the superconducting transition temperature \Tc \cite{DuanNature2021,RaymondJPSJ2021}.
    So \UTe might be close to the case of UPt$_3$, where neutron scattering studies also detected only antiferromagnetic correlations \cite{Aeppli1988}: nevertheless, some  theories attempted to explain the observed $E_{2u}$ superconducting ground state of UPt$_3$ with a pairing mechanism based on pure antiferromagnetic fluctuations \cite{NomotoPRL2016}.

    In \UTe, theoretical models of the spin-triplet superconducting state have explored three main scenarios, starting from band structure calculations:
    i) ferromagnetic fluctuations, winning over antiferromagnetic fluctuations for specific values of the normal state parameters (exchange constants, Coulomb repulsion...)\cite{Xu2019, Ishizuka2021}; ii) local (intra unit-cell) ferromagnetic correlations between the nearest neighbour uranium ions \cite{Shishidou2021, HazraArXiv2022}; iii) pure antiferromagnetic fluctuations, thanks to multi-orbital degrees of freedom \cite{ChenArXiv2021} or peculiar conditions on the $\bm{Q}$-dependent susceptibility and the underlying Fermi-surface \cite{KreiselPRB2022}.

    Our detailed thermodynamic measurements demonstrate that the field reinforcement of \Hc in UTe$_2$ arises from a new superconducting phase, driven by a pairing mechanism different from that controlling the low field superconducting phase.
    Conversely, we also found that in the low field phase, pairing is suppressed by magnetic fields applied along the easy magnetisation axis.
    All together, we show that the superconducting symmetry may change from a spin-triplet to a spin-singlet state as a function of magnetic field on approaching the metamagnetic instability. All of these new features are solid experimental inputs challenging the theoretical scenarios for the superconducting pairing in \UTe.

    The present paper is organised as follows:
    in the next section~\ref{sec:overview}, we give a rapid overview of the main experimental results.
    Section~\ref{sec:experiment} gives the experimental details.
    The results are presented in section~\ref{sec:results} and analysed in section~\ref{sec:analysis}.
    Section~\ref{sec:discussion} discuss their impact and the open questions.
    More data and details on the methods used for the analysis can be found in the Supplemental Material \cite{Supplemental}.

\section{Overview of the main results}
\label{sec:overview}
    This work reveals a thermodynamic phase transition between the low and high-field reinforced superconducting phases, for magnetic fields applied along the hard magnetisation \baxis.
    This has been made possible thanks to the use of two bulk thermodynamic probes, namely specific heat up to 36~T and thermal dilatation/magnetostriction measurements up to 30~T.
    We also uncover a drastic alteration of the specific heat anomaly along \Hc between the two phases, an additional support to the change of pairing mechanism between the two superconducting phases.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/FigOV_1_PD_Hb.pdf}
	\caption{Top: Bulk superconducting phase diagram of \UTe for a magnetic field $H$ along the \baxis.
	The specific heat measurements reveal a phase transition line between a low-field (\lftransi) (yellow) and high-field (\hftransi) (green) superconducting state.
	As a function of magnetic field, superconductivity is suppressed at 34.75~T by a first-order metamagnetic transition marking the entrance in a partly polarised magnetic phase (magenta).
	Bottom: Specific heat divided by temperature as a function of temperature at different magnetic fields: (left) At 12~T a sharp anomaly marks the superconducting transition to the \lftransi phase. (middle) On cooling a broad hump-like transition occurs at the transition to the \hftransi phase above a sharp low temperature transition to the \lftransi superconducting phase. (right) At 26~T only the hump-like transition to the \hftransi superconducting phase is observed.
	}
	\label{fig:SimplePhaseDiagram}
\end{figure}

    Fig.~\ref{fig:SimplePhaseDiagram} is a summary of this main result, showing the phase diagram for field along the \baxis, with the two superconducting phases (labelled \lftransi and \hftransi for "low-field and "high-field"), and the specific heat anomalies at the different phase transitions.
    Phase transitions between superconducting phases of different symmetries are rare, however such a change in the specific heat anomaly between the different phases is, to the best of our knowledge, unique.
    Analysis is required to reveal that it suggests different spin-states characterising these two phases, with a \hftransi phase in strong interplay with the metamagnetic transition occurring, for field along the \baxis, at \Hm$=34.75$~T.
    The result which is rather counter-intuitive is that the \lftransi phase would be spin-triplet, and the \hftransi phase spin-singlet, most likely triggered by the development of antiferromagnetic correlations on approaching \Hm.
    Hence, this would be a direct consequence of the competing pairing interactions (or of the change of dominant finite $\bm{Q}$-vector in the magnetic excitation spectrum) predicted to occur in \UTe \cite{Xu2019, Ishizuka2021, KreiselPRB2022}.

\begin{figure}
	\centering
	\includegraphics[width=.9\linewidth]{PicturesMain/FigOV_2_Hc2_Ha.pdf}
	\caption{Top: \Hc along the \aaxis. Dashed line : orbital limit deduced from the measured initial slope at \Tc.
	Bottom: zoom on the very low field behaviour, close to Tc, showing the very strong negative curvature responsible for the deviation of \Hc along $a$ from the orbital limitation.
    This anomalous behaviour signs a strong suppression of the pairing strength for fields applied along the easy magnetisation axis, an effect also observed in the ferromagnetic superconductors.}
	\label{fig:thePbWithHc2_along_a}
\end{figure}

    Another surprise uncovered by these specific heat measurements is that \HcT appears to be anomalous not only along the \baxis, with the field-reinforced \hftransi phase, but also along the $c$ and most importantly along the easy \aaxis.
    This had been completely overlooked up to now, but determination of \HcT by specific heat reveals a very strong negative curvature of \Hc along the \aaxis very close to \Tc, and an initial slope four times larger than initially thought.
    Hence, a new mechanism is required to explain the very strong deviation of \HcT from the linear behaviour expected below \Tc/2 (see Fig.~\ref{fig:thePbWithHc2_along_a}):
    we show that a strong paramagnetic limitation, anyway excluded by recent NMR Knight-shift measurements \cite{FujibayashiJPSJ2022}, would not be sufficient to explain this singular temperature dependence of \Hc.
    Instead, it points to a strong suppression of the pairing strength along the easy axis.
    In the case of \UTe, it could arise from at least two different sources.
    Suppression of (hypothetical) ferromagnetic fluctuations by fields along the easy axis would lead to such a decrease of the pairing strength, a  mechanism similar to that in ferromagnetic superconductors \cite{Hattori2012, Wu2017}, but also working for paramagnetic systems \cite{Mineev2017}.
    Or the strong sensitivity of finite-$\bm{Q}$ (spin-triplet) pairing \cite{ChenArXiv2021, KreiselPRB2022} to Fermi-surface instabilities, already revealed in \UTe at rather low magnetic fields along the \aaxis \cite{NiuPRL2020}, might also lead to pairing strength suppression.

    Hence, both results along the hard \baxis and the easy \aaxis bring new elements on the possible symmetry states in \UTe and on the competing pairing mechanisms.
    It enlightens the stunning superconducting properties of \UTe, and uncovers key features which should guide future theory developments.
    Indeed, understanding the mechanisms leading to the strong field dependence of the pairing strength in \UTe is an issue barely touched by current theoretical models, but certainly central for this system and clearly of major interest for the whole field of unconventional superconductivity.

\section{Experimental Details}
\label{sec:experiment}
\subsection{Single crystal growth and samples}

	Different single crystals of \UTe from three different batches have been studied by specific heat measurements and magnetostriction/thermal expansion measurements.
	All single crystals were prepared by the chemical vapour transport method with iodine as transport medium.
	A starting ratio of U:Te $=$ 2:3 has been used, and the quartz ampules have been heated slowly up to a final temperature gradient of 1060$^\circ$C/1000~$^\circ$C, which was kept for 18 days.
	The ampules have been slowly cooled down to ambient temperature during 70~hours.
 	The first sample ($\#$1) has a superconducting transition temperature \Tc of 1.45 K, with a mass of 12.3~mg. This sample is from the same batch than those studied in Refs.~\cite{Knebel2019, NiuPRR2020}. It has been studied by specific heat and magnetostriction.  Samples, \diese2 and \diese3 with masses of 5.6 mg and 27 $\mu$g respectively, are from a different batch with a critical temperature around 1.85~K.  In addition, thermal expansion has been measured on a fourth crystal (sample \diese4). Samples  \diese2, \diese3 and  \diese4 show almost similar values of the specific heat in the normal state, and of the specific heat jumps at $T_c \approx 1.85$~K as well as of the residual value $\gamma_0$. (We found $\gamma_0 \approx 0.03$~J K$^{-2}$mol$^{-1}$ at 0.1~K, or 0.011~J K$^{-2}$mol$^{-1}$ extrapolated from above 0.3~K.)

\subsection{Specific heat measurements}
	The specific heat of two samples ($\#$1 and $\#$2) has been measured by a quasi-adiabatic relaxation method in a dilution refrigerator up to 15~T in a superconducting magnet and down to 100~mK. Small heat pulses of maximum 1$\%$ of the temperature $T$ (0.5$\%$ in the superconducting transition) were applied to the sample.
	The specific heat $C$ is extracted from the temperature response of the sample during the whole pulse sequence.
	Down to the lowest temperatures, only one relaxation time was measured in the exponential decay.
	The addenda have been measured separately.
	They represent $8\%$ of the total measured specific heat at 2~K and $2\%$ at 100~mK.
	Mainly temperature sweeps were performed, but also some field sweeps for the transitions at the lowest temperatures (between 100~mK and 200~mK).

    To align the samples in the magnetic field, we used a piezoelectric rotator allowing a rotation over 90$^{\circ}$ in a plane parallel to the field, and a goniometer allowing a $\pm$3$^{\circ}$ rotation perpendicular to the plane.
	Furthermore, the set-up is rigid, so that the torque between magnetic field and the anisotropic magnetisation of the sample could not induce a misalignment.\\

	 Sample (\diese3) of 27 $\mu$g coming from the same batch as sample \diese2, has been measured with an ac specific heat technique in a $^3$He refrigerator down to 600~mK, and up to 36~T in the M9 magnet at the high magnetic field laboratory LNCMI in Grenoble.
	 Details of the specific heat set-up are shown in the Supplemental Materials of Ref.~\onlinecite{MichonNature2019}.
	 For fields up to 18.5~T, the ac calorimetry has been performed using a 20~T superconducting magnet (M2) in combination with a $^3$He refrigerator down to 400~mK.
	 An attocube piezorotator allowed for a rotation in the ($b$,$c$) plane.\\

	 The value of the critical temperature \Tc is extracted from the specific heat transition by a fit to an ideal jump broadened by a Gaussian distribution of critical temperatures: \Tc corresponds to the centre of this distribution.
     This model, which reproduced the data very well, allows to extract directly other parameters like the width and jump of the specific heat at the transition (more details in the Supplemental Material Section V.B.).

\subsection{Linear magnetostriction and thermal expansion measurements}
    The linear magnetostriction $\Delta L_b/L_b$ of \UTe has been measured on sample \diese1.
    In addition, we measured the linear thermal expansion at constant magnetic field on a second single crystal \diese4.

    These measurements have been performed using a high resolution capacitance dilatometer \cite{Kuechler2012}.
	The capacitance has been determined using an Andeen Hagerling capacitance bridge AH2550A.
	High magnetic field experiments have been performed using the 30~T magnet M10 of the high magnetic field laboratory LNCMI Grenoble.
	Due to the limited diameter of the magnet it has been only possible to measure the length change $\Delta L_b$ of the sample parallel to the magnetic field applied along the $b$ axis of the crystal.
    The magnetic field has been swept with a maximal rate of 100 G/sec to avoid eddy currents heating.
	The dilatometer was positioned at the end of a silver cold finger in a $^3$He cryostat, with a base temperature near 370~mK. A RuO$_2$ thermometer and a heater were fixed directly on the dilatometer. Temperature sweeps at fixed magnetic field have been performed with maximal heating rates of 0.1~K/min.

	Additional thermal expansion measurements have been performed in a superconducting magnet up to 13~T using a dilution refrigerator in CEA Grenoble.

\section{Results}
\label{sec:results}
	\subsection{Magnetic field \texorpdfstring{$H \parallel b$\xspace}{H parallel b}}
    \subsubsection{Specific heat}
	In zero field, all samples studied exhibit a single sharp superconducting transition (width $\Delta T_c \approx 20$ to $38$~mK), with a large jump at the superconducting transition (up to $\Delta$\Cp/\Cp$\approx$~1.85), emphasising the high quality and homogeneity of the samples.
	Figure~\ref{fig:RawCpLF} shows $C/T$ as a function of temperature measured on sample \diese2 for several magnetic fields $H \parallel b$\xspace up to 15~T.
	In zero field the transition is extremely sharp at $T_c = 1.85$~K.
	Under magnetic field it remains sharp with a pronounced jump $\Delta C/T_c$ up to 18.5~T (see Fig.~\ref{fig:RawCpLF_HF}-a), so the transitions are easily followed under field.
	The specific heat measurements on the first \UTe samples displayed an upturn and a large residual term of $C/T$ at low temperatures \cite{Metz2019a}.
	Both became smaller with improved sample quality.
	Our measurements on crystal \diese2 show indeed only a small residual term, and an upturn shifted to lower temperatures compared to samples with lower \Tc (see also in Supplemental Material \cite{Supplemental} Section V.A.).
	This agrees with more recent work claiming that residual term and upturn are extrinsic to \UTe \cite{Cairns2020,AokiJPCM2022,Rosa2022}.

\begin{figure}[pt]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig3_RawCpLF.pdf}
	\caption{Temperature dependence of the specific heat $C/T$ at different magnetic fields $H \parallel b$ from 0~T to 15~T measured on sample \diese2. }
	\label{fig:RawCpLF}
\end{figure}

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig4_RawCpLF_HF.pdf}
	\caption{(a) Temperature dependence of $C/T$ measured on sample \diese3 for fields $H \parallel b$\xspace from 17 to 18.5~T.  A second wide transition appears above the sharp low temperature transition. (b) $C/T$ for the different magnetic fields up to 35.5~T, which is above the metamagnetic transition. }
	\label{fig:RawCpLF_HF}
\end{figure}

	Remarkably, above 15~T for $H \parallel b$\xspace, $C/T$ shows two transitions.
    A second wide transition appears for fields above 15~T (350~mK width at 18~T), which is well detached above $H \gtrsim 17$~T from the sharp transition observed at lower fields (Fig.~\ref{fig:RawCpLF_HF}-a).
	We follow this second wide \hftransi transition (Fig.~\ref{fig:RawCpLF_HF}-b) up to the metamagnetic transition
	\cite{Knafo2019, MiyakeA2019, RanNatPhys2019}: its width and height remain roughly constant with field.

    This is best seen through a Gaussian analysis of the temperature dependence of \CT, allowing to deconvolute broadening effects.
    The specific heat jump at the transition from the \lftransi to \hftransi phase displays a marked drop compared to its amplitude for the transition between the normal and the superconducting \lftransi phase (see Fig.~\ref{fig:Cp_Jump_Width}).
    Essentially, it goes down to the same level as that of the wide transition of the \hftransi phase, which remains roughly constant up to \Hm .
    Hence, as expected, the emergence of the \hftransi transition goes along with a redistribution of entropy between both phases.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig5_Cp_Jump_Width.pdf}
	\caption{(a)Jump of the transitions $\Delta C/T_c$ as a function of field, for $H \parallel b$\xspace, determined on sample \diese2 and \diese3 for the \lftransi and \hftransi transitions.
	(b) Width $\Delta T$ of the transitions as a function of fields. $\Delta T$ is equal to 2.35 times the standard deviation of the Gaussian distribution of T$_c$ used to fit the transition.}
	\label{fig:Cp_Jump_Width}
\end{figure}

	These high-field measurements shown on Fig.~\ref{fig:RawCpLF_HF} were performed on  crystal \diese3, which is from the same batch as \diese2 in Fig.~\ref{fig:RawCpLF}.
	As observed already from resistivity measurements in the field-reinforced superconducting phase, the \Tc of this broad anomaly is increasing with field, except very close to \Hm where the transition temperature decreases slightly.
	This may be due to a slight misalignment of the sample in the high-field experiments (no sample rotation), or it could be intrinsic.
	Above \Hm, the broad anomaly abruptly disappears.
	This \hftransi transition is most likely the expected bulk signature of the field-reinforced superconducting phase observed in transport properties for the same field direction \cite{Knebel2019, RanNatPhys2019}.
	We will see in the next section that we have other experimental evidence that the sample displays bulk superconductivity in this high-field region.

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig6_RawCpFieldSweepLF_HF.pdf}
	\caption{$C/T$ versus field for $H \parallel b$\xspace on sample \diese3 at different temperatures between 0.7~K and 1.86~K. The \lftransi superconducting transition and the peak at \Hm$=34.75$ T are clearly visible, but the \hftransi superconducting transition reported on Fig.~\ref{fig:RawCpLF_HF} appears only as a very broad anomaly.}
	\label{fig:RawCpFieldSweepLF_HF}
\end{figure}

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig7_RawCpHm.pdf}
	\caption{$C/T$ versus field ($H \parallel b$\xspace) at the metamagnetic transition for sample \diese3 at different temperatures between 0.7~K and 1.86~K. The hysteresis is independent of the field sweep rate. On the left panel, the arrows indicate the direction of the field sweeps.}
	\label{fig:RawCpHm}
\end{figure}

	Figure \ref{fig:RawCpFieldSweepLF_HF} displays the field dependence along the $b$ axis of $C/T$ up to 36~T.
	On approaching the first order metamagnetic transition at \Hm,
	\CT shows a strong increase, with a large drop (of order 25 \%) above \Hm, and a hysteresis (see Fig.~\ref{fig:RawCpHm}), independent of the sweep rate.
	The drop of the specific heat is sharpest at the lowest temperature, with a width of 0.25~T.
	However, a possible interplay between the superconducting and metamagnetic transition at this temperature may influence the shape of the anomaly.
	The width of the hysteresis decreases linearly with increasing temperature, starting from 0.17~T at 0.7~K.
    This behaviour of \CT at \Hm agrees qualitatively with previous measurements \cite{Imajo2019, MiyakeA2021} performed in pulsed magnetic fields (see comparison in Supplemental Material \cite{Supplemental} Section IV).
    As regards the superconducting transitions below \Hm, the sharp \lftransi transition is well observed in the field sweeps (see Fig.~\ref{fig:RawCpFieldSweepLF_HF}).
    However, the \hftransi transition observed in temperature scans appear as a very broad anomaly, and it is noticeable only by comparison with curves at different temperatures.
    This arises due to the combination of an already large \Tc distribution at fixed field, and an almost vertical \Hc, so that this \hftransi transition appears extremely broad as a function of field (see Supplemental Material \cite{Supplemental} Section V.D.).

\subsubsection{Complete phase diagram}
	Since samples \diese2 and \diese3 come from the same batch and have essentially the same \Tc at 0~T (1.847~K and 1.845~K respectively), we use the specific heat measurements on these two samples to establish the complete superconducting phase diagrams for all field directions from 0 to 36 T.
	It is shown on Fig.~\ref{fig:PhaseDiagCp}.
	As underlined already in section \ref{sec:overview}, the most important result is that for $H \parallel b$\xspace, two superconducting phases are clearly present, with a point around $H=15$~T and $T=1$~K where the three transition lines should join.

	The low field (\lftransi) superconducting phase corresponds to the sharp \lftransi transition that can be followed from 0~T up to 18.5~T.
	The high field (\hftransi) phase corresponds to the broad \hftransi transition emerging above 15~T, and followed up to \Hm (Figures~\ref{fig:RawCpLF}-\ref{fig:RawCpLF_HF}).
	For magnetic fields slightly above 15~T, the two transitions overlap, thus it is difficult to determine unambiguously \Tc for the \hftransi phase transition. Knowing that the jump and width of these transitions seem to have a negligible field-dependence (Fig.~\ref{fig:Cp_Jump_Width}), we fixed them to extract $T_c(H)$  between 15 and 17~T (empty crosses on Fig.~\ref{fig:PhaseDiagCp}). (More details on the analysis and the way we determine the transition temperatures are given in the Supplemental Material \cite{Supplemental} Section V.B.)
	Hence these values need to be taken with some caution, and there is yet no experimental determination of the precise nature of the crossing of the \Hc lines of the \lftransi and \hftransi phases: they could merge with a sharp change of slope, or could be tangential.
	Moreover, we looked thoroughly but found no sign of a fourth transition line inside the \lftransi phase, either on temperature or field sweeps, as would be expected for a crossing of second order phase transitions \cite{YipPRB1991, Braithwaite2019}.
	This absence of a fourth line however could also arise from a further broadening of the \hftransi transition inside the \lftransi phase (for field below 15~T).

	Nevertheless, the phase diagram clearly demonstrates the existence of at least two different superconducting phases for $H \parallel b$\xspace:
	comparison with the resistivity measurements, and the observation of vortex pinning in the (\hftransi) phase by linear magnetostriction ( see next paragraph) show that it is also a superconducting phase.

\begin{figure}[pt]
	\centering
	\includegraphics[width=0.9\linewidth]{PicturesMain/Fig8_PhaseDiagCp.pdf}
	\caption{Phase diagram up to 36 T for $H$ applied along the three crystallographic axis, established with specific measurements on samples \diese2 (below 15~T) and \diese3 (above 15~T).
	Blue circles: $H \parallel a$\xspace, green triangles: $H \parallel c$\xspace, red symbols: $H \parallel b$\xspace. Red squares are \Tc from the sharp \lftransi transitions, crosses are \Tc from the broad \hftransi transitions.
	For empty crosses, \Tc could only be determined by fixing the width and jump of the \hftransi transition. Magenta hexagons: \Hm determined on sample \diese3 by specific heat.}
	\label{fig:PhaseDiagCp}
\end{figure}

\subsubsection{Linear Magnetostriction}
	The longitudinal linear magnetostriction $\Delta L_b(H)/L_b$ along the \baxis for $H \parallel b$\xspace of sample \diese1 is shown on Fig.~\ref{fig:RawdL}.
	In the normal state ($T =2$~K), the linear magnetostriction is negative and shows roughly a $H^2$ field-dependence, as usually observed in paramagnetic metals (Fig.~\ref{fig:RawdL}(a)).
	This is in agreement with the low field measurements of Ref.~\cite{ThomasPRB2021} and with the very recent measurements in pulsed magnetic fields, which show a strong negative jump of the linear magnetostriction at the metamagnetic transition \cite{MiyakeA2022}.
	Following Maxwell's relations, the negative sign of $\Delta L_b/L_b$  indicates that under uniaxial stress applied along the \baxis, the susceptibility $\chi_b$ along this axis should increase.
	We note that such an increase of $\chi_b$ has already been observed under hydrostatic pressure \cite{Li2021}.

	\begin{figure}
		\centering
		\includegraphics[width=1\linewidth]{PicturesMain/Fig9_RawdL.pdf}
		\caption{(a) Longitudinal linear magnetostriction $\Delta L_b/L_b$ of sample \diese1 along the $b$ axis for a field applied along the \baxis, in the superconducting state at 0.35~K and in the paramagnetic state at 2~K. (b) Linear magnetostiction in the superconducting state at different temperatures after subtraction of the paramagnetic contribution measured at 2~K. The grey arrows indicate the direction of the field sweep. Black vertical arrows mark the closing and reopening of the hysteresis as a function of field, which corresponds to the irreversibility field $H_{\rm irr}$ of the vortex motion in the superconductor.}
		\label{fig:RawdL}
	\end{figure}

	In the superconducting state at $T = 0.35$~K, the linear magnetostriction shows a very pronounced hysteretic behaviour.
	In Fig.~\ref{fig:RawdL}(b) we display the additional contribution to the linear magnetostriction $\Delta L_b^s$ which appears in the superconducting state.
	It is obtained from the measured linear magnetostriction at fixed temperature in the superconducting state after subtraction of the paramagnetic contribution measured at 2~K.
	The linear magnetostriction in the superconducting state is very large and shows a strong hysteresis with a fishtail-like behaviour both below $\approx 15$~T and above $\approx 20$~T.

	This irreversible magnetostriction appears very similar to the behaviour of the magnetisation in the mixed state of type~II superconductors with strong vortex pinning.

	In  the critical state model \cite{Bean1964}, magnetic flux penetration or expulsion, when increasing (or decreasing) the field, is impeded by vortex pinning.
	This imposes a field gradient at the sample surface, perpendicular to the applied field, controlled by the critical current density.
	Hence, if magnetic flux lines are trapped by the action of pinning forces, equal but opposite forces will act on the lattice.
	Thus, the length change of the crystal $\Delta L_b/L_b$ is proportional to $\frac{H\Delta M}{E}C_\nu$, where $\Delta M$ is the non-equilibrium part of the magnetisation, $E$ the Young modulus, and $C_\nu$ a constant depending on the Poisson's ratio \cite{Eremenko1999}.

	In \UTe, at the lowest temperature, the hysteresis in the linear magnetostriction vanishes above $H_{\rm irr,1} \approx 14.5$~T, and it opens again above $H_{irr,2} \approx 21$~T, being maximal at 30~T, which is the highest field we could reach in the experiment.
	On increasing the temperature, the lower field $H_{\rm irr,1}$ decreases, while the upper field $H_{\rm irr,2}$ increases and above 0.6~K, it exceeds the achievable field range.

  	\begin{figure}[pt]
	\centering
	\includegraphics[width=0.9\linewidth]{PicturesMain/Fig10_PhaseDiadL.pdf}
	\caption{Phase diagram magnetic field vs temperature for field $H \parallel b$ for sample \diese1 with a \Tc of 1.45~K.
	Blue circles: irreversibility fields determined from the magnetostriction measurements (see Fig.~\ref{fig:RawdL}); green diamonds: \Hc from specific heat measurements performed on the same sample; red triangles: \Hc from thermal conductivity on a sample from the same batch (see Supplemental Material \cite{Supplemental} Section I), and (dashed line) \Hc determined from the resistivity $\rho (H)$ of another sample from the same batch
	(data from Ref.~\cite{Knebel2019}).
	Magenta line: \Hm determined from resistivity (circles, \cite{Knebel2019}).

	}
	\label{fig:PhaseDiagdL}
	\end{figure}

    In Fig.~\ref{fig:PhaseDiagdL} we display the irreversibility field $H_{\rm irr}$ of sample \diese1, ($T_c \approx 1.5$~K) determined from the magnetostriction measurements.
    We compare this field with \Hc determined from specific heat measurements (only below 15~T on this sample), with thermal conductivity measurements (see section I in the Supplemental Material \cite{Supplemental}) and also with resistivity measurements on a sample from the same batch with similar \Tc (published in \cite{Knebel2019}).
    In the \lftransi phase, thermodynamic as well as thermal and electrical measurements are in good agreement, with only small quantitative differences on the amplitude of the negative curvature of \HcT in the 0-5~T field range.
    Obviously, flux pinning in the sample at low magnetic fields is rather strong, as found by low field-low temperature magnetisation measurements \cite{Paulsen2021}, and the irreversibility  field follows the upper critical field: $H_{irr,1} \sim H_{c2}$.
    In the field-reinforced \hftransi superconducting phase above 15~T, the difference between $H_{irr,2}$ and $H_{c2}$ gets more pronounced: there is a broader reversible regime, between \Hc and $H_{irr,2}$.
    Actually, the observation of the irreversible magnetostriction due to the flux pinning between 20 and 30~T is a first proof that the \hftransi phase delimited by the broad specific heat anomaly is indeed a bulk, field-reinforced superconducting phase.

    In the Supplemental Material (\cite{Supplemental} Section II)  we also show the linear thermal expansion $\Delta L_b(T)/L_b$ as a function of temperature measured on sample \diese4, which has a \Tc = 1.82~K similar to that of samples \diese2 and \diese3 studied by specific heat.
    On sample \diese4, we could observe both \Hc, which shows up as a kink in the temperature dependence of the sample length (hence a jump of the thermal expansion $\Delta L_b/L_b$), as well as the entrance of the irreversible regime of the magnetisation, controlled by vortex pinning.
    In the \hftransi phase, as opposed to the \lftransi phase, (\cite{Supplemental} Section II), the  temperatures for \Hc and for the emergence of the irreversible regime are well separated.  This suggests a decrease of the pinning strength in the \hftransi phase.

\subsection{\texorpdfstring{\Hc}{Hc2} close to \texorpdfstring{\Tc}{Tc}, along all directions}

    The upper critical field \Hc in \UTe present anomalies not only due to the presence of these two superconducting phases for $H \parallel b$\xspace:
    it also has an anomalous temperature dependence close to \Tc for all directions.
    The initial slope of \Hc at \Tc (\dHc)  has its largest value for $H \parallel b$\xspace.
    From our specific heat measurements, we determine \dHc$\approx-34$ T/K, which is larger than the values obtained from electrical transport measurements.
    Generally speaking, in electrical transport measurements, filamentary superconductivity easily leads to slightly higher \Tc values.
    However this filamentary superconductivity is fragile under small fields, hence it is usual to observe a small positive curvature of \Hc determined by resistivity.
    In \UTe, it may mask the negative curvature of \Hc clearly visible in Fig.~\ref{fig:PhaseDiagCp} close to \Tc.

    Figure~\ref{fig:PhaseDiagCp} also shows the upper critical field \Hc of sample \diese2 along the $a$ and $c$ axes.
    Similar to the previous reports \cite{RanScience2019, Knebel2019}, \Hc is strongly anisotropic and extrapolates to 9~T for the $a$ axis, 15~T for the $c$ axis.
	Along the \caxis, \Hc is almost linear on a large temperature range, and has the lowest slope near \Tc, \dHc$\approx-7.5$ T/K.

\begin{figure}[pt]
	\centering
	\includegraphics[width=1.\linewidth]{PicturesMain/Fig11_Hc2LF.pdf}
	\caption{Upper critical field \Hc as a function of temperature determined from specific heat measurements on sample \diese2 for the $a$, $b$ and $c$ axes at very low fields, close to Tc. The dashed line is a linear fit following the initial slope of \Hc along the easy magnetisation \aaxis.
	It is close to the value along the \baxis, however followed by a very strong negative curvature showing up at very low fields (0.4~T is about $\mu_0 H_{c2}(0)/20$ for $H \parallel a$\xspace).}
	\label{fig:Hc2LF}
\end{figure}

    The most striking feature appears along the \aaxis: it seems also linear and with a slope at \Tc comparable to that along the \caxis, when examined on this broad field and temperature scale.
    A closer look very near \Tc yields a contrasted view, showing that \Hc along this easy axis has a large initial slope, followed by a strong negative curvature.
    Figure~\ref{fig:Hc2LF} displays a zoom on the very low field region ($H<0.4~T$):
    it can be seen that the initial slope \dHc is of order -20~T/K for $H \parallel a$\xspace, much larger than the values determined by resistivity measurements.

    Actually, this effect was probably already present in Ref.~\onlinecite{Kittaka2020} reporting also specific heat measurements, but was barely discussed.
    We have confirmed this anomalous behaviour (the very large slope followed by a very strong curvature), on three different samples measured also by specific heat (see Supplemental Material \cite{Supplemental} Section V.C.).
    Discussion of the possible origin of this very strong curvature along the \aaxis is of course central for the question of the symmetry state of the superconducting order parameter in \UTe and of the pairing mechanism.

	\subsection{Normal phase specific heat}
	It is now well established, notably from the very large specific heat jump at \Tc, that \UTe is in a strong-coupling regime.
	The most natural explanation for the field-reinforced superconducting phase is that the superconducting pairing itself is enhanced under fields along the \baxis \cite{RanScience2019, Knebel2019}.
	As discussed already for the ferromagnetic superconductors, in the strong coupling regime, such a field dependence should be reflected also in the normal state Sommerfeld specific heat coefficient $\gamma$ \cite{Wu2017, MineevAnnPhys2020}.
	In \UTe, it is not so easy to determine the Sommerfeld coefficient, because \CT remains strongly temperature dependent almost down to \Tc, and there is no simple way to analyse this temperature dependence:
	this arises mainly from a marked anomaly in the specific heat with a maximum at a temperature $T^{\ast} \approx 12$~K \cite{Willa2021a}
	(see Supplemental Material \cite{Supplemental} Section III).
	Hence, as already reported in the literature \cite{MiyakeA2019, Imajo2019, MiyakeA2021}, we show in Fig.~\ref{fig:Cp_T_FieldSweep} $C(H)/T$ normalised to $C/T$ in zero field, at $T=1.86$~K for fields applied along the three crystallographic directions.
	As the temperature dependence of \CT becomes weak at this temperature, it can be considered as a reasonable estimation of the $\gamma$ value, at least as long as the applied field is not too close from \Hm.
	Indeed, a specific heat anomaly at 12~K in zero field, and attributed to magnetic fluctuations \cite{Willa2021a}, is shifted to much lower temperatures on approaching \Hm for $H \parallel b$\xspace \cite{Willa2021a,AokiJPCM2022}.
	Hence, the jump of \CT at \Hm is not necessarily a jump of the Sommerfeld coefficient, it can also have a large contribution from the suppression of thermal magnetic fluctuations.

	\begin{figure}[pt]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig12_Cp_T_FieldSweep.pdf}
	\caption{$C/T$ normalised by its value at zero field, as a function of field for $H$ along $a$,$b$ and $c$ axis at 1.86 K. Measurements below 15 T have been done on sample \diese1, by temperature sweeps. Measurements for $H \parallel b$\xspace above 15~T have been done on sample \diese3 by a field sweep at 1.86 K. Insert shows a zoom on the measurements below 15 T, lines are guides to the eyes.}
	\label{fig:Cp_T_FieldSweep}
\end{figure}

	Along the \baxis, $\frac{C}{T}(H)$ increases  under field, and even more strongly on approaching \Hm.
	It is similar to already published data \cite{MiyakeA2019, Imajo2019}, with some quantitative differences notably above \Hm, where the ac technique in static fields probably allows for more precision.
	By contrast, for the field along the $a$ or $c$ axis, $C/T(H)$ has a much more complex behaviour.
	It is known that along the \aaxis, Lifshitz anomalies appear around 5 T and possibly 9~T \cite{NiuPRL2020}.
	They are indeed visible on our measurements (see Supplemental Material \cite{Supplemental} Section III).
	In addition, we observe anomalies also as pronounced maxima along the $a$ and $c$ axes, respectively at 0.5 T and 1.5 T, whose origin is still unclear.
	Hence, until better understood, it is difficult to rely on the field dependence of $C/T$ to discuss the behaviour of the pairing strength with field.

\section{Analysis}
\label{sec:analysis}
	\subsection{\texorpdfstring{\Hc}{Hc2} in the  \texorpdfstring{\lftransi}{LF} phase}

    Some discussion is first required on the behaviour of \Hc close to \Tc.
    Indeed, the observed negative curvature of \Hc along the $b$ and $a$ axes might seem to contradict the common belief that \UTe is a \pwave superconductor with a \dvector perpendicular to the easy \aaxis.

    Indeed, the very strong curvature of \Hc along the \aaxis may suggest a severe paramagnetic limitation.
    This is at odds with the value of $H_{c2}(0) \sim 9$~T, which is much larger than the weak-coupling paramagnetic limit of around 3.5~T (for a gyromagnetic factor $g=2$ with \Tc=1.85~K).
    In effect, the negative curvature is so concentrated close to \Tc that it requires a large value of $g\sim 3.2$ (in the weak-coupling limit) to match the initial deviation from linearity of \Hc along the \aaxis, leading to a saturation of \Hco at 2.25~T at low temperatures (see also Fig.~S13 of the Supplemental Material \cite{Supplemental} Section V.F.).
    In other words, \HcT along the easy axis does not follow at all the temperature dependence of an upper critical field solely controlled by paramagnetic and orbital limitations \cite{WerthamerPR1966_3}.
    Adding strong-coupling effects will not change anything to this problem, except for requiring even larger values of $g$ to reproduce the initial negative curvature (\cite{Supplemental} Section V.F.).

    Before discussing further this anomalous temperature dependence of \Hc along the \aaxis, let us examine the new larger value of the initial slope derived from these bulk measurements (see Fig.~\ref{fig:PhaseDiagCp}).
    Previous measurements of the lower critical field \Hcl \cite{Paulsen2021} had already pointed out an inconsistency between \Hcl and \Hc along the \aaxis.
    From the value of \Hcl, it was deduced (\cite{Paulsen2021}-Supplemental Material) that \dHc should be -16~T/K (instead of -6.6~T/K deduced from resistivity measurements, for a sample with \Tc$\sim 1.45$~K), which rescales up to -20.4~T for a \Tc of 1.85~K.
    Hence, the observed value of -20~T/K is now fully consistent with \Hcl along the \aaxis.

    As regards the same comparison of \Hcl with \Hc along the \caxis, the agreement is also very good.
    Again, after rescaling \Tc, from the values of \dHcl along the \caxis in \cite{Paulsen2021}, we expect a value of \dHc of -7.6~T/K, whereas our experiments yield to -7.5~T/K.

    This has to be contrasted with the case for $H \parallel b$\xspace.
    As stated already in \cite{Paulsen2021}, the anisotropy of \Hcl between the $b$ and $c$ axes at \Tc is very small.
    Hence \dHc should be roughly equal (and of the order of -8~T/K) in both directions, whereas \dHc is largest along the \baxis and of the order of -34~T/K (see Table~S1 in Supplemental Material \cite{Supplemental} Section V.E.). In addition, the temperature dependence of \Hc along the \caxis and the \baxis also cannot be reproduced by any combination of paramagnetic and orbital limitations:
    for $H \parallel c$\xspace, \HcT is too linear with even a small positive curvature close to \Tc, and for $H \parallel b$\xspace, the situation is similar to that along the \aaxis, however with smaller deviations visible mainly above 8~T (Fig.~S13 in the Supplemental Material \cite{Supplemental} Section V.F.).

    The only way we can imagine to understand these anomalies regarding the value of the slope at \Tc (for $H \parallel b$\xspace) and the temperature dependence of \Hc along all directions, is a field dependent pairing strength.
    This has been observed in ferromagnetic superconductors \cite{Wu2017} and been already proposed for \UTe \cite{RanScience2019, Knebel2019}.
    If we call $\lambda$ the strong-coupling parameter controlling this pairing strength, it has to be field-dependent in all directions.
    We can rely on \Hcl, where the effects of such a field dependence have negligible effects \cite{Paulsen2021} to fix the average Fermi velocities controlling \dHc for $H \parallel b$\xspace, without the contribution of the field-dependent pairing strength.
    Along the $a$ and $c$ axes, the agreement between \dHcl and \dHc shows that $\frac{d\lambda(H)}{dH} \sim 0$ in zero field (at \Tc).

    For the estimation of $\lambda(H)$ in the different directions, there are general constraints which are model-independent.
    First of all, along the \baxis, the discrepancy between \dHcl and \dHc can only be reconciled with an increase of the pairing strength:
    increasing \dHc requires $\frac{d\lambda(H)}{dH}(H=0) > 0$.
    Hence, we expect an increase $\lambda(H)$ in the (\lftransi) phase.

    For $H \parallel c$\xspace, starting with $\frac{d\lambda(H)}{dH}(0) = 0$, the small positive curvature and the very linear behaviour also requires an increase of $\lambda(H)$, whatever the strong-coupling model and the mechanisms (orbital and/or not paramagnetic limitation).

    The situation is not so straightforward for $H \parallel a$\xspace.
    Indeed, the deviation from linearity observed very close to \Tc can arise most naturally by a strong suppression of $\lambda(H)$ with an \Hc otherwise purely orbitally-limited.
    But in principle, it could also be explained with an increasing $\lambda(H)$ (see Supplemental section ).
    Quantitatively, this second scenario seems unrealistic, requiring a very strong paramagnetic limitation completely at odds with the NMR Knight-shift measurements \cite{FujibayashiJPSJ2022}.

     Conversely, the most natural scenario of a suppressed pairing strength with \Hc controlled mainly by the orbital limit is consistent with the Knight-shift measurements , yielding essentially no change of the Knight-shift along the \aaxis, and only minute changes along the $b$ and $c$ axes.
    It should be stressed here that NMR measurements are performed at fixed field.
    They extract the change of electronic spin susceptibility (when there is a mean to estimate the "reference level" in the normal state) across \Tc from the temperature variation of the Knight-shift.
    Hence, this measurement is not directly influenced by the field dependence of the pairing strength, as opposed to considerations on the violation of the paramagnetic limit on \Hc.

    In the following, we will concentrate on this second hypothesis: spin-triplet superconductivity with an equal-spin pairing state (ESP state) characterised by a \dvector perpendicular to the \aaxis.
    It is more reasonable as regards values of the gyromagnetic factor, and it is well supported by other experiments.
    Following the NMR results \cite{FujibayashiJPSJ2022}, we also assumed a negligible paramagnetic limitation in the (\lftransi) phase along the $b$ and $c$ axes.
    We used the same strong-coupling model as in \cite{Wu2017,Knebel2019}, which mimics \Hc for a spin-triplet superconductor with a calculation for an \swave superconductor with $g=0$.
    Average Fermi velocities have been adjusted to match the slopes of \dHc calculated from \dHcl and the thermodynamic critical field (Table~S1 in the Supplemental Material \cite{Supplemental} Section V.E.).

    The result for $\lambda(H)$ in the three crystallographic directions is reported in Fig.~\ref{fig:Lambda_H}.
    The increase of $\lambda$ along the $c$ and even $b$ axes is modest in the (\lftransi) phase, at most 10\%, whereas a strong suppression (factor 2 between \Tc and $T \rightarrow 0$) is required for $H \parallel a$\xspace.
    The inset of Fig~\ref{fig:Lambda_H} shows that despite the strong suppression of the pairing for $H \parallel a$\xspace, $\frac{d\lambda}{dH}=0$ in zero field.
    As explained above, this originates from the absence of correction on the initial slope of \Hc$\parallel a$ at \Tc, imposed in the fitting procedure.
    Conversely in \ucoge, the slope of \Hc along the easy axis is strongly suppressed already at \Tc.
    This difference between \UTe and \ucoge is fully consistent with predictions for paramagnetic and ferromagnetic superconductors respectively, where $\frac{d\lambda}{dH}$ due to the suppression of "ferromagnetic" fluctuations is proportional to $M_z\frac{\partial M_z}{\partial H}$ ($M_z$ being the magnetisation along the easy axis) \cite{Mineev2017}.

    The different parameters of the fit are summarised in Table~S2 of the Supplemental Material \cite{Supplemental} Section V.F.
    It is to be noted that the lowest Fermi velocities (\vF) are along the $b$ and $c$ axes and highest along \aaxis.
    This matches qualitatively the anisotropy found in transport measurements between the different axes \cite{Eo2021}.
    Quantitatively, \Hc depends on an average of the Fermi velocities perpendicular to the applied field direction, weighted by the pairing strength.

    As a consequence for example, we never succeeded to compare quantitatively anisotropies of \Hc in UPt$_3$ or in URu$_2$Si$_2$ with the detailed determination of their respective Fermi surfaces by quantum oscillations \cite{McMullan_2008, Bastien2019}, even though the order of magnitude of the orbital limitation is consistent with measured effective Fermi velocities.
    This is probably even more acute if subtle $\bm{Q}$-dependent paring is responsible the specific pairing state realised in the system, as could well be the case in \UTe \cite{Xu2019,Ishizuka2021,KreiselPRB2022,ChenArXiv2021}.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig13_Lambda_H.pdf}
	\caption{$\lambda(H)$ determined from \Hc measured on sample \diese2. $\lambda$ was set to 1 at 0~T, $g$ to 0 for the \lftransi phase. Red squares are for $H \parallel b$\xspace, green triangles for $H \parallel c$\xspace and blue circles for $H \parallel a$\xspace. For the \hftransi phase, red crosses symbolise \llambdaH determined with $g=2$, and grey crosses with $g=0$. A zoom on fields below 1~T is shown in the insert.}
	\label{fig:Lambda_H}
\end{figure}

	\subsection{\texorpdfstring{\Hc}{Hc2} in the  \texorpdfstring{\hftransi}{HF} phase}

    The field-reinforced \hftransi phase is a central property of \UTe and from the very start, two mechanisms known to lead to a reinforced \Hc in conventional superconductors were evoked and rejected \cite{RanScience2019}.
    A Lebed mechanism \cite{Lebed2020}, could lead to a similar behaviour of \Hc, but it would not explain the double transition revealed in this work, and quantitatively, it does not fit with the overall 3D character of the system \cite{AokiJPCM2022}.
    A Jaccarino-Peter mechanism \cite{Jaccarino1962} associated to a competing singlet (or at least, strongly paramagnetically limited phase) could lead to such a field-induced phase \cite{Meul1984}, but it would require a very large orbital limit, hence necessarily a competing singlet phase with a larger \Tc in zero field than the present \lftransi phase.

    So, theoretical models have proposed that it could arise from a symmetry change of the order parameter \cite{IshizukaPRL2019, Shishidou2021}.
    The main idea is that for a spin-triplet superconducting state arising from ferromagnetic fluctuations along the easy magnetisation axis, at low fields, the \dvector should be perpendicular to this $a$ axis. Hence a $B_{3u}$ (or more generally $B_{3u}+iB_{1u}$ state under field \cite{AokiJPCM2022}) is favoured at low fields.
    By contrast, for strong fields along the \baxis, a rotation of the \dvector is expected toward a $B_{2u}$ state (or $B_{2u}+iA_{u}$), to minimise the component of the \dvector along the \baxis.
    Such a symmetry change would imply a phase transition somewhere between the low and high field regimes, which had not been detected until the present specific heat measurements.
    However, this change of \dvector orientation alone could only explain that an initial paramagnetic limitation is exceeded thanks to the new orientation of the \dvector. It will not explain the positive \dHcT observed in the \hftransi phase between 15 and 30~T.

    Conversely, empirical explanations of the reinforcement of \Hc focused on field-induced enhancement of the pairing (positive $\frac{d\lambda}{dH}$)  \cite{RanScience2019, Knebel2019}, which is compatible with, but does not require an additional phase transition.
    Indeed, even rather sharp upturns observed on \Hc extracted from electrical transport can be reproduced with a smooth continuous increase of $\lambda(H)$ \cite{AokiJPCM2022}.

    This is not the case anymore with the present experiments, showing not only that there is a field-induced thermodynamic phase transition between two different superconducting states, but also that the pairing mechanisms driving these phases is likely different.
    This is first seen from the phase diagram of Fig.~\ref{fig:PhaseDiagCp}, where the \lftransi phase is only weakly field-enhanced, extrapolating to 20~T at zero temperature, whereas the \hftransi phase requires a much steeper field increase:
    there is more than just a symmetry change and rotation of the \dvector for the transition of the \lftransi to \hftransi phases.
    The change of the superconducting anomaly is another support to the change of pairing mechanism.
    It differs strongly from the case of UPt$_3$ \cite {HasselbachPRL1989} or more recently of CeRh$_2$As$_2$ \cite{KhimScience2021, Landaeta2022}, where no change is observed on the shape of the specific heat transition along \Hc when switching from the low to the high field superconducting phases.
    In \UTe, the specific heat anomaly for both phases is markedly different, with a very broad anomaly in the \hftransi phase whereas that of the \lftransi phase remains remarkably sharp (see Fig.~\ref{fig:RawCpLF_HF}).

    Most likely, pairing is reinforced at high fields due to the closer proximity to \Hm.
    The phase diagram for $H \parallel b$\xspace under pressure \cite{Knebel2020} shows that the \hftransi phase could well be the same as the pressure-induced higher \Tc superconducting phase \cite{Braithwaite2019,RanPRB2020,Aoki2020,ThomasPRB2021}.
    There are two main theoretical proposals for this "high-\Tc" pressure-induced phase.
    The first \cite{Shishidou2021} is that it is a $B_{2u}$ phase, having no component of the \dvector along the \baxis, with a pairing mechanism controlled by local ferromagnetic correlations.
    The second is that it is a spin-singlet ($A_g$) phase \cite{Ishizuka2021}, induced by antiferromagnetic correlations becoming dominant over ferromagnetic fluctuations under pressure.

    The first proposal is the most natural one, explaining the different phases with a single pairing mechanism and transitions imposed by the Zeeman coupling of the \dvector with the applied field.
    In this framework, nothing should be changed to the procedure for the evaluation of the field dependent pairing strength $\lambda(H)$ between the \lftransi and \hftransi phases.
    The result is reported also on Fig.~\ref{fig:Lambda_H} (grey crosses), displaying a cusp at $\sim 15$~T, but a rather weak increase (30\% between zero field and \Hm) of the pairing strength, far from the factor$\approx 2$ observed on \CT at 1.8~K (Fig.~\ref{fig:Cp_T_FieldSweep}).

\begin{figure}[pb]
	\centering
	\includegraphics[width=1\linewidth]{PicturesMain/Fig14_FitCpHF.pdf}
	\caption{The \hftransi transition at 24 T, $H \parallel b$\xspace, measured by specific heat.
	Lines are the transitions calculated for the spin-singlet (g=2) or spin triplet (g=1) superconducting state in the \hftransi phase, and for the measured distribution of \Hm (1$\sigma$) or a distribution twice as large (2$\sigma$).
	For the spin-triplet state, even with a distribution of \Hm twice larger than given by our measurement, we fail to reproduce the broadening.
	By contrast, for the same larger distribution of \Hm, the agreement is good for the spin-singlet scenario.}
	\label{fig:FitCpHF}
\end{figure}

    Using the other proposal of a spin-singlet superconducting order parameter for the \hftransi phase seems paradoxical, but yields interesting results.
    Once again, with a field dependent pairing, the weak-coupling paramagnetic limit is easily exceeded thanks to the effective increase of the "zero-field" \Tc, and strong coupling effects.
    Under pressure, this singlet phase can explain the strong paramagnetic limitation observed along the \aaxis \cite{Knebel2019,Aoki2020a}.
    Here, it leads to a stronger field dependence of the pairing strength (red crosses on Fig.~\ref{fig:Lambda_H}), required to overcome the saturation of \Hc (at fixed $\lambda$, see Fig.~S14 in the Supplemental Material \cite{Supplemental} Section V.F.).
    Estimation of $\lambda(H)$ has been done for a gyromagnetic factor of 2, and (arbitrarily) with the same energy scale $\Omega$ than for the \lftransi phase, considering that both mechanisms should have similar characteristics in order to lead to similar critical temperatures.
    Using different values of $\Omega$ (but the same $g$-factor) changes little to the following analysis.

    Indeed, this spin-singlet scenario can explain a large part of the broadening observed for the superconducting transition in the \hftransi phase.
    Considering that the field dependence of $\lambda$ might come from the proximity to \Hm, a simple hypothesis is that $\lambda$ is a function of $\frac{H}{H_m}$.
    Then, a dispersion of \Hm controlling the broadened specific heat anomaly reported on Fig.~\ref{fig:RawCpHm} for the metamagnetic transition, will translate into a distribution of \Tc, hence to a new mechanism for the broadening of the superconducting transition.
    From the calculation of \Hc at fixed $\lambda$ used to extract the field dependence of the pairing, we can as well determine $T_c = \varphi\left(H, \lambda\left( \frac{H}{H_m}\right)\right)$.
    This allows to determine the effect of the distribution of \Hm (deduced from an analysis of the specific heat at \Hm), on the specific heat anomalies of the \hftransi phase, according to the two different determinations of $\lambda(H)$ (see Fig.~\ref{fig:FitCpHF}).

 	\begin{figure}[pt]
 		\centering
 		\includegraphics[width=1\linewidth]{PicturesMain/Fig15_FitCp_theta.pdf}
 		\caption{Temperature dependence of $C/T$ measured on sample \diese3 $H \parallel b$\xspace at 18.5~T for different angles in the ($b$,$c$) plane. Dashed lines are the transitions calculated from a distribution of \Hm, controlled by its angular dependence \cite{RanNatPhys2019} and a finite mosaicity of 3$^{\circ}$ in the sample.}
 		\label{fig:FitCp_theta}
 	\end{figure}

    With this hypothesis, the dispersion of \Hm explains half the width of the superconducting transition in the spin-singlet scenario (it fits well the observed anomaly by doubling the \Hm dispersion) but it fails in the spin-triplet case.
    A simple analysis (see Supplemental Material \cite{Supplemental} Section V.G.) shows that the key difference between the two scenarios is the much larger value of $\frac{\partial T_c}{\partial H}\vert_{\lambda}$ imposed by the paramagnetic limitation.

    Within this scheme, the dispersion of \Hm is found to have a negligible influence on the \lftransi transition, so that it does give a hint for why the two superconducting phases could be marked by such different specific heat anomalies.

    Moreover, when the sample is misaligned in the ($b$,$c$) plane (see Fig.~\ref{fig:FitCp_theta} for an applied field of 18.5~T, and more data in the Supplemental Material \cite{Supplemental} Section V.H), the transition shifts to lower temperatures and the amplitude of the jump decreases.
    At 15$^{\circ}$ the \hftransi transition almost disappears.
    Taking only into account  angular dependence of \Hm \cite{RanNatPhys2019}, and a hypothetical mosaicity of 3$^{\circ}$ in our crystal, we can roughly reproduce the huge broadening of the anomaly at finite angles, with the same dependence of \Tc on \Hm (dashed dotted lines on Fig.~\ref{fig:FitCp_theta}), and otherwise a constant ideal specific heat jump.
    This is another support for this explanation of the large broadening of the specific heat anomaly, relying on the sensitivity of the pairing strength to $\frac{H}{H_m}$, and a paramagnetic limitation of the upper critical field in the \hftransi superconducting phase.

\section{Perspectives}
\label{sec:discussion}
	A clear result from this work is the requirement of a field-dependent pairing strength along all directions of the applied field, as shown by the anomalous temperature dependence of \Hc along the $a$, $b$ and $c$ axes already close to \Tc.
	The (most likely) strong decrease of the pairing strength along the \aaxis is reminiscent of the results on \ucoge along its easy magnetisation axis, and at first sight, it seems best compatible with a pairing mechanism involving true ferromagnetic fluctuations.
	For paramagnetic systems like \UTe, this mechanism also leads to a vanishing slope of $\frac{d \lambda}{d H}$ along the easy axis, consistent with the good agreement found between \Hcl measurements and the present determination of \dHc.

	There are also several theoretical studies showing that spin-triplet pairing is possible with finite momentum magnetic fluctuations \cite{KreiselPRB2022}, or with only local ferromagnetic correlations within a unit cell \cite{Shishidou2021}.
	The field dependence of such mechanisms hasn't been explored.
	However, the Fermi surface instability observed at 6~T along the easy axis \cite{NiuPRL2020} could play a role if $\bm{Q}$-dependent pairing is important.
	Hence, the origin of this strong suppression of the pairing along the easy axis in \UTe requires deeper and more detailed  investigations, before strong conclusions can be drawn.
	A detailed characterisation of the Fermi surface of \UTe is notably a mandatory step.
	In any case, hasty comparison  with ferromagnetic superconductors could be misleading.

	This is even clearer regarding the results along the hard \baxis.
	From this work, we can now claim the existence of two different bulk superconducting phases in \UTe already at ambient pressure, when a magnetic field is applied along the \baxis above 15~T.
	Analysis of these results has been done looking for the simplest possible explanations.
	In the \lftransi phase, we confirm that contradiction between \Hcl and \Hc values close to \Tc requires a steep increase of the pairing strength for $H \parallel b$\xspace already at low fields ($\frac{d \lambda}{d H} > 0$ at \Tc).
	In the \hftransi phase, we have found support for a paramagnetically limited \Hc along the \baxis, as it could explain a large part of the strong broadening of the specific heat anomaly in the \hftransi phase, and the still increased broadening when turning away from the \baxis in the ($b$,$c$) plane.

	In the course of revising this paper, we became aware of a new work on high-field NMR, recovering a similar phase diagram as reported here, but identifying the \hftransi phase as a spin-triplet $A_u+iB_{2u}$ state \cite{KinjoArXiv2022}.
	This arises from Knight-shift measurements in the \hftransi phase, showing no detectable changes across \Tc.
	We note that due to the field dependence of the pairing strength, these measurements in the \hftransi phase are all performed at values of $\frac{H}{H_{c2}^{eff}(0)}$ close to 1, where $H_{c2}^{eff}(0)$ is the effective value of \Hco, corresponding to the value of the pairing strength ($\lambda(H)$ ) at the field ($H$) of the measurement (see Fig.~S14 of the supplement \cite{Supplemental}).
	At these large field values (with respect to $H_{c2}^{eff}(0)$), there is probably little change to expect for the Knight shift,  whatever the spin-state.

	Naturally, many open questions still remain.
	A first one concerns the large difference between the irreversibility line observed on the magnetostriction and the specific heat anomaly:
	it suggests a strong increase of the reversibility region in this phase, like the fact that the resistive transition only goes to zero when the bulk transition ends.
	This is similar to previous observations in \ucoge \cite{Wu2018}, and  we still have no explanation for this phenomenon.

	Another important point concerns the order of the different transition lines and the precise slopes of the lines at the tricritical point.
	Indeed, as for CeRh$_2$As$_2$ \cite{KhimScience2021}, in the case of a spin-triplet to spin-singlet transition, a first order transition is expected.
	In our specific heat measurements, we did not detect hysteresis effects.
	The only features visible on Fig.~\ref{fig:Cp_Jump_Width} are a smaller jump of \CT for the \lftransi to \hftransi transition than along \Hc in the \lftransi phase, as well as a slightly smaller width.
	This slight narrowing leaves open the possibility that the transition from \lftransi to \hftransi phases could be weakly first order.
	There are many cases in condensed matter physics, where first order transitions lead to negligible hysteresis:  see e.g. Ref.~\onlinecite{GuillouNatCom2018}, or the old example of the $^3$He melting curve which can therefore be used for the transmission of the Provisional Low Temperature Scale (PLTS-2000) between 0.9~mK and 1~K.
	Nevertheless, in \UTe this point requires further experimental investigations.

    If this transition is first order, of course, the question of the tricritical point is solved.
    If it is not, it remains an issue to determine if there is an additional transition line within the \lftransi superconducting phase, and whether or not the three transition lines determined in this work join with different slopes, or if the \Hc line has no change of slope (only a very strong positive curvature) at the tricritical point.

    Note that the entrance into the \hftransi phase along \Hc cannot be done in a mixed singlet-triplet superconducting phase \cite{YanasePrivate}:
    it would require, like for the chiral superconducting state  \cite{Hayes2021}, a double transition which is not observed.

   The interplay between ferromagnetic and antiferromagnetic fluctuations \cite{Xu2019, Ishizuka2021} could lead to competing pairing interactions \cite{Ishizuka2021}.
    This competition could be central both for the pressure and the field-induced phases of \UTe.
   At ambient pressure, at the opposite of CeRh$_2$As$_2$, it would lead to a paradoxical spin-singlet phase at high fields, most likely driven by strong antiferromagnetic correlations on approaching the metamagnetic transition.
   Under pressure, this high-field phase would become the highest \Tc phase with the lowering of the metamagnetic field along the \baxis, whereas the pure spin-triplet phase would survive essentially for large enough fields along the easy \aaxis.
   \UTe would then be probably the first system where two competing pairing mechanism of similar strength exist, and can be arbitrarily tuned by field or pressure.
   Finding in the same system both spin states, controlled "at will" by magnetic field or pressure is an ideal case to challenge theoretical models and understand which conditions allow for the emergence of spin-triplet superconductivity

\section{Acknowledgements}
We thank Y. Yanase, M. Houzet, M. Zhitomirsky and V. Mineev, for very fruitful discussions, and  A. Miyake and S. Imajo for kindly providing their data points.
We got financial support from the CEA Exploratory program TOPOHALL, the French National Agency for Research ANR within the project FRESCO No. ANR-20-CE30-0020 and FETTOM ANR-19-CE30-0037,
and from the JSPS programs KAKENHI (JP19H00646, JP20K20889, JP20H00130, JP20KK0061, JP20K03854, JP22H04933).
We acknowledge support of the LNCMI-CNRS, member the European
Magnetic Field Laboratory (EMFL), and from the Laboratoire d\'excellence LANEF (ANR-10-LABX-51-01).

\section{Thermal conductivity \texorpdfstring{$H \parallel b$\xspace}{H parallel b}}
\label{secS:thermalconductivity}
Thermal conductivity $\kappa$ measurements have been performed on a sample from the same batch as sample \diese1.
The temperature sweeps have been measured on a home-made dilution refrigerator with a base temperature of 100~mK and a superconducting magnet with field up to 16~T using a standard "one heater-two thermometers" setup. The temperature dependence of $\kappa/T$ for different magnetic fields up to 16~T is represented in Fig.~ \ref{figS:KappaUTe2}. $\kappa/T$ shows a broad maximum at around 3~K. At low field, there is a clear increase in $\kappa/T$ just below \Tc, which is suppressed by increasing the field. Such an increase below \Tc has also been observed in other systems such as CeCoIn$_5$ or YBCO, and attributed to a suppression of the inelastic scattering of heat carriers (electrons and phonons, respectively) by the opening of the superconducting gap. In the case of \UTe, the enhanced conductivity below \Tc is likely due to an increase of the electronic mean free path due to the opening of the superconducting gap. At higher field ($\mu_0 H>3$~T), entrance in the superconducting state is marked by a rapid decrease of the thermal conductivity, usually attributed to Andreev scattering on the vortex cores.

\begin{figure}[hb]
\includegraphics[width=0.4\textwidth]{PicturesSupp/FigS1_KappaUTe2.pdf}
\caption{Temperature dependence of $\kappa/T$ in \UTe with $H \parallel b$\xspace between 0 and  16~T (every Tesla). Traces have been shifted for clarity. The superconducting temperature \Tc is represented by vertical arrows.}
\label{figS:KappaUTe2}
\end{figure}

\clearpage
\section{Thermal expansion for field \texorpdfstring{$H \parallel b$\xspace}{H parallel b}}
\label{secS:thermalexpansion}

In addition to the linear magnetostriction measurements discussed in the main text, we performed longitudinal thermal expansion  measurements (length change in the field direction) on sample \diese4 as a function of temperature at fixed magnetic fields up to 29.5~T along the \baxis in the high magnetic field laboratory LNCMI and in addition using a superconducting magnet up to 13~T in the Pheliqs laboratory.

Figures \ref{figS:RawDilatation_PhaseDiag} (a) and (b) show the temperature dependence of the relative length change along the $b$ axis $\Delta L_b/L_b$ for different magnetic fields.
The data displayed in panel (a) are obtained by cooling from the normal state to the lowest temperature in a field 2~T below the target final field, then increasing the field at the lowest temperature up to the final field, heating at that field above $T_c$ and cooling again.
Increasing the field at low temperature induces a non-equilibrium magnetisation inside the sample due to the flux pinning, while an equilibrium flux distribution occurs in the final-field cooled sweep.
$H_{\rm irr}$ marks the onset of the irreversible magnetisation regime.

The data in panel (b) are measured in a superconducting magnet using a dilution refrigerator.
They are obtained from a similar procedure however starting by zero field cooling, then ramping the field up to its target value, heating up to the normal state and measuring the field-cooled length change. Again, the irreversibility line can be clearly determined.

In panel (c) we show the temperature dependence of the upper critical field \Hc and of the irreversibility field $H_{\rm irr}$ determined from the cycles described above. Data for $H \leq 13$~T are from the experiments performed in a superconducting magnet, while for $H > 13$~T they stem from the high-field experiments in LNCMI. While in the LF superconducting phase $H_{\rm irr}$ is lower but very close to  $H_{\rm c2}$, in the HF superconducting phase $H_{\rm irr}$ is far lower in temperature than \Hc, indicating a wide reversible region (with low pinning ) behind \Hc in this state.

\begin{figure}[hb]
\includegraphics[width=0.8\textwidth]{PicturesSupp/FigS2_RawDilatation_PhaseDiag.pdf}
\caption{Temperature dependence of the linear thermal expansion $\Delta L_b/L_b$ in \UTe with $H \parallel b$\xspace at different magnetic fields measured in (a) the LNCMI Grenoble, and (b) using a superconducting magnet. The arrows indicate the direction of the temperature sweep, see text for details. The temperature of the irreversibility field $H_{\rm irr}$ and the upper critical field \Hc are indicated by arrows and vertical bars, respectively.
Panel (c) shows the phase diagram obtained from the thermal expansion measurements. While at low field ($H < 15$~T) $H_{\rm irr} \approx H_{\rm c2}$, at high fields the irreversibility line is much lower in temperature. }
\label{figS:RawDilatation_PhaseDiag}
\end{figure}

\clearpage
\section{specific heat: normal phase}

    In general, the leading term of the specific heat of a heavy-fermion system at low temperatures ($T \ll T_F$, with $T_F$ being the effective Fermi temperature), far from any quantum criticality, is linear in temperature: $C \propto \gamma T$.
    The Sommerfeld coefficient $\gamma$  is proportional to the density of states at the Fermi level, which is, strongly renormalised compared to the free electron gas.
    In UTe$_2$, at low temperature, the competition between different natures of magnetic fluctuations (ferromagnetic or antiferromagnetic) as well as the role of valence fluctuations due to the interplay between U$^{3+}$ and U$^{4+}$ configurations may occur.
   The situation is even more complex in this system, because electronic correlations play the unusual role of driving the system from an insulating toward a metallic state \cite{Aoki2019, Xu2019, Shick2019, Ishizuka2019}.

    As a consequence, even close  to \Tc,  $C/T$ is not the sum $\gamma + \beta T^2$ with $\beta T^2$  the phonon contribution far below the Debye temperature. An additional  contribution is observed, likely coupled to the "Schottky-like" anomaly detected at $T^\star \approx 12$~K \cite{Willa2021a}.
    In Fig.~\ref{figS:ct_vs_T2}(a,b) shows $C/T$ vs $T^2$ for sample \diese1 and sample \diese2.
    The phonon contribution has been calculated following the Debye model with a Debye temperature of $\theta_D = 180$~K \cite{Willa2021a}.
    Clearly, the phonon contribution is low compared to the measured specific heat and cannot reproduce the strong temperature increase of \CT at zero magnetic field.
    Under magnetic field, the anomaly at $T^\star$ shifts to higher temperatures for $H \parallel a$\xspace and lower temperatures for $H \parallel b$\xspace.

   Accordingly, \CT strongly depends on the magnitude of the applied field and on its direction.
   Hence the low temperature Sommerfeld coefficient $\gamma$ cannot be determined properly from a direct analysis of the temperature dependence of \CT.
   We tried, as done previously \cite{MiyakeA2019,Imajo2019}, to follow the field evolution of \CT at the lowest possible temperature, so as to be as close as possible to the value of $\gamma$.

    At 15~T for $H \parallel a$\xspace, $C/T$ decreases with temperature (down to 0.3~K), which is quite unusual (Fig.~\ref{figS:ct_vs_T2}-c).
    Thermoelectric power measurements have revealed the presence of several Lifshitz transitions in this field direction \cite{Niu2020}.
    In the main text, the insert of Figure~\ref{fig:Cp_T_FieldSweep} shows that the field dependence of $C/T$ at 1.8~K has a minimum followed by a maximum close to 5~T and 9~T respectively.

 \begin{figure}[b]
	\centering
	\includegraphics[scale=.6]{PicturesSupp/FigS3_CT_vs_T2.pdf}
 	\caption{Temperature dependence of $C/T$ (sample \diese1) for fields applied along the three crystallographic directions.
 	(a) $C/T$ as a function of $T^2$, measured on sample \diese1.
 	No linear behaviour is seen.
 	At 15~T for $H \parallel a$\xspace, the temperature dependence is drastically suppressed compared to measurement at 0~T.
 	At 15~T for $H \parallel b$\xspace, the temperature dependence is slightly larger compared to zero field.
 	The doted line corresponds to the phonon contribution calculated from estimation of the Debye temperature by high-temperature measurements on other samples \cite{Willa2021a}, in addition to a constant Fermi-liquid contribution).
 	(b) Same data for $H \parallel b$\xspace on sample \diese2.
 	The anomalous magnetic contribution seems more pronounced than for sample \diese1.
 	(c) $C/T$ at low temperatures at 15 T along the three axis measured on sample \diese1. The superconducting transition at $\sim 0.5$~K is visible for $H \parallel b$\xspace.}
	\label{figS:ct_vs_T2}

\end{figure}

    Moreover, field sweeps have been performed up to 31~T for $H \parallel a$\xspace , on a sample with a \Tc of 1.45~K coming from the same batch as sample \diese1. Measurements of \CT are displayed in Fig.~\ref{figS:CT_vs_H_para_a}.
    A minimum (around 6 T) followed by a maximum (around 8 T) are visible for fields above the superconducting transition.
    At higher fields, a change of slope occurs for fields between 17~T and 22~T, depending on the temperature.
    We can follow these three anomalies in addition to the superconducting transition, and establish the phase diagram shown in Fig.~\ref{figS:PhaseDiagLifschitz}.
    The temperature dependence and the order of magnitude of the field where they occur are similar to those obtained from thermoelectric power measurements. The lower transition (at around 5-6~T) was clearly identified as a Lifshitz transition, the origin of the two others is less clear \cite{Niu2020}.

    Regarding the present specific heat data, the origin of the pronounced maxima of \CT observed at $\sim 1$~T for $H \parallel a$\xspace and $\sim 1.5$~T for $H \parallel c$\xspace (Fig.~\ref{fig:Cp_T_FieldSweep}) is also not identified.

\begin{figure}[ht]
   \begin{minipage}{0.48\textwidth}
    \centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS4_CT_vs_H_para_a.pdf}
	\caption{$C/T$ measured for field sweeps at different temperatures with $H \parallel a$\xspace.
	Measurements done on a sample coming from the same batch as sample \diese1, with a \Tc of 1.5 K.}
	\label{figS:CT_vs_H_para_a}
   \end{minipage}\hfill
   \begin{minipage}{0.48\textwidth}
     \centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS5_PhaseDiagLifschitz.pdf}
	\caption{Phase diagram $H \parallel a$\xspace up to 31 T. Black triangles represent the superconducting transitions, blue circles are the minima of $C/T(H)$.
	Red squares represent the maxima of $C/T(H)$.
	Green pentagons represent the inflexion point of $C/T$ observed on field sweeps.
	The dashed lines are the corresponding transitions measured by thermoelectric power in ref~\cite{Niu2020}.}
	\label{figS:PhaseDiagLifschitz}
   \end{minipage}
\end{figure}

\clearpage
\section{Specific heat: metamagnetic transition}
\label{secS:SuppCpHm}
    As regards the metamagnetic transition, we could measure precisely the field dependence of the anomaly, and check that it did not depend on field sweep rate, varied between $\pm 350$ and $\pm 50$~Gauss/sec.
    At this sweep rate, we also did not detect any magnetocaloric effect.
    Hence, the present continuous field sweep measurement show unambiguously that there is a jump at \Hm, slightly broadened by a distribution of \Hm. At 0.7~K, \CT decreases 55~mJK$^{-2}$mol$^{-1}$ and the width of the transition is 0.25~T.
    This distribution of \Hm possibly comes from the strong sensitivity of \Hm to pressure \cite{Knebel2019,MiyakeA2022} and most likely to stress (crystal defects could then generate this \Hm distribution).

    Figure~\ref{figS:CompCT_Imajo} shows the comparison of $C/T - C/T(H=0)$ near $T=1.8$~K determined from our experiment performed on sample \diese3 and experiments performed in pulsed field.
    Ref.~\onlinecite{Imajo2019} reports specific heat experiments in highly stabilized fields, using the long pulsed fields facility at ISSP.
    In Ref.~\onlinecite{MiyakeA2021} the Sommerfeld coefficient $\gamma$ has been determined from the magnetization measurements $M(T)$ under pulsed fields, using Maxwell's relation for $H \ne H_m$ as
    $(\partial \gamma/\partial H)_T = (\partial ^2M/\partial T^2)_H$, and  using the Clausius-Clapeyron relation for the first order transition:
    $\mu_0 dH_m/dT = -\Delta S/\Delta M$ \cite{MiyakeA2021} to get the jump
    $\Delta\gamma = \Delta S/T$ at \Hm.
    The last analysis indicated a discontinuous jump of  $\Delta \gamma = -30$~mJK$^{-2}$mol$^{-1}$ at $H_M$ for $H \parallel b$\xspace, which is lower than that obtained in the present experiment. However, despite some quantitative differences (e.g. the absolute variation of $C/T - C/T(H=0)$ is larger in both pulsed field experiments) the general behaviour is similar: an increase with $H$ when approaching \Hm and a drop at \Hm followed by a strong decrease. A similar field dependence has been observed for the $A$ coefficient of the electrical resistivity, albeit without a clear jump above  \Hm\cite{Knafo2019a, Knafo2021a}.

\begin{figure}[ht]
   \begin{minipage}{0.48\textwidth}
        \centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS6_CompCT_Imajo.pdf}
\caption{Cyan triangles: specific heat measurements done in pulsed fields in ref~\cite{Imajo2019}. Orange circles: $\gamma(H)-\gamma(H=0)$ determined from the magnetisation measurements through thermodynamic relations in ref~\cite{MiyakeA2021}.
Red line: our measurements.}
\label{figS:CompCT_Imajo}
   \end{minipage}\hfill
   \begin{minipage}{0.48\textwidth}
     \centering
	\includegraphics[scale=0.45]{PicturesSupp/FigS7_param_transi_Hm.pdf}
	\caption{(a) Jump (drop) of $C/T$ at \Hm, (b) \Hm and (c) width at half-height ($2.35 \sigma$) of the transitions as a function of the temperature, for the up and down sweeps.}
	\label{figS:param_transi_Hm}
   \end{minipage}
\end{figure}

To extract the position and the width of the specific heat at the metamagnetic transition, it has been fitted, like a broadened second-order phase transition, with the same Gaussian model as the superconducting transition.
Fields replace temperatures, and \Hm replaces \Tc. The parameters extracted are displayed in Figure~\ref{figS:param_transi_Hm}.
\Hm for the up sweeps are roughly constant at 34.75 T, and for the down sweeps \Hm increases with the temperature. This goes along with the trend toward a closing of the hysteresis cycles with increasing temperatures, in agreement with the observations from resistivity measurements \cite{Knafo2019a,Niu2020b}.
The width of the transitions increases abruptly from 0.24 T to 0.45 T between 0.7 K and 0.97 K and then stays constant with temperature.
The jump of \CT at \Hm strongly decreases on cooling from  0.97 K to 0.7 K, otherwise, above 0.97K,  it decreases on warming.
This anomaly at 0.7~K might be due to the presence of the \hftransi transition, which is wide enough in field and temperature to influence the drop of \CT at \Hm.

\clearpage
\section{Specific heat: superconducting phase}

\subsection{Measurements in zero field}
All measurements of \CT in \UTe display an upturn below 0.1K, and an extrapolated (from temperatures above the upturn) residual term at $T=0$ which was quite large in the first measurements \cite{RanScience2019,Aoki2019a,Metz2019}.
More recent studies are claiming that the residual term and the upturn are extrinsic to \UTe \cite{Cairns2020,AokiReview2022}.
Our measurements on different samples in Fig.~\ref{figS:CompCpSamplesBT}(a-b) show diverse behaviours at low temperatures.
The upturn is not monotonously correlated to \Tc, however, it is strongly reduced on our best samples.
The residual term seems to be more systematically decreasing with the \Tc increase, well in the trend reported in \cite{AokiReview2022}.
In any case, these measurements do agree with an extrinsic nature of these anomalies.
Note also that on sample \diese2, the entropy balance is perfectly satisfied at \Tc, within experimental errors (better than 1\%).
The low temperature upturn plays a negligible role (see Fig.~\ref{figS:CompCpSamplesBT}(c)).

\begin{figure}[ht]
   \begin{minipage}{0.48\textwidth}
        \centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS8_CompCpSamplesBT.pdf}
	(a)
    \end{minipage}\hfill
   \begin{minipage}{0.48\textwidth}
    \centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS8_CompCpSamplesBT_loglog.pdf}
  (b)
    \end{minipage}
    \vspace{0pt}
   \begin{minipage}{0.48\textwidth}
	\includegraphics[scale=0.5]{PicturesSupp/FigS8_entropy.pdf}
   (c)
    \end{minipage}
 	\caption{$C/T$ as a function of temperature at zero field and low temperatures for different samples.
	Samples \diese1 and \diese2 are the one presented in this article. Sample neutron is a large sample of 241~mg used to perform the neutron diffusion experiments in ref~\cite{Knafo2021,RaymondJPSJ2021}.
	(a) linear scale, (b) logarithmic scale, emphasising the strong decrease of the residual term when the critical superconducting temperature in enhanced.
	(c) Entropy calculated for samples \diese1 and \diese2, showing the bad balance for \diese1 and the very good one for \diese2.}
	\label{figS:CompCpSamplesBT}
\end{figure}
\clearpage
\subsection{Gaussian model for the specific heat anomaly}
\label{secS:Gaussian}
A simple hypothesis is that broadening of the specific heat transition is controlled by a Gaussian $T_c$ distribution of the form:
\begin{equation}
    p(T_c) = \frac{1}{\sigma \sqrt{2 \pi}} \exp \left( -\frac{1}{2}\left( \frac{T_c - T_{c0}}{\sigma} \right)^2\right)
\end{equation}

For the specific heat, or any additive quantity we can then write that:
\begin{equation}
	C/T = \int_{-\infty}^{\infty} p(T_c)  C/T(T,T_c)  dT_c
	\label{equS:CpDistrib}
\end{equation}

The simplest expression for $C/T(T,T_c)$ is a constant $\gamma$ term above \Tc, a jump at \Tc followed by a constant positive slope below \Tc. If both the slope and the jump are independent of \Tc, this amounts to:
\begin{equation}
	C/T(T,T_c)  = \gamma + \theta(T_c - T) \left( \deltaCT + \alpha (T-T_c) \right)
	\label{equS:CpIdeal}
\end{equation}

Hence, for the total specific heat:
\begin{equation}
	C/T(T)  = \gamma + (\deltaCT + \alpha (T-T_{c0}))\left[\frac{1}{2}  - \frac{1}{2 } \text{erf} \left( \frac{T-T_{c0}}{\sigma} \right) \right]  - \alpha \frac{\sigma} {\sqrt{2 \pi}} \exp \left(-\frac{1}{2}\left( \frac{T - T_{c0}}{\sigma}  \right)^2 \right)
	\label{equS:CpAnalytique}
\end{equation}

This is model is fine for the zero field transition, where $\deltaCT$ is independent of $T_c$.
However, under magnetic field, the broadening of the transition may correspond to a distribution of slopes of \Hc (proportional to $T_c$ for clean type II superconductors).
We can expect that $\deltaCT$ will be suppressed with field, with a decrease controlled by $H/H_{c2}(0)$.
Hence $\deltaCT$ will not be constant within the broadened transition.
More simply, we can assume that the jump will be suppressed like $T_c(H) / T_c(0)$.
The problem is therefore to relate $T_c(H)$ and $T_c(0)$, or more precisely, to get the $T_c(0)$ corresponding to a given $T_c(H)$.
Then we could take for a model of the transition that $\deltaCT$ is proportional, within the transition, to $T_c(H) / T_c(0)$.
A simple way to find this relation is to assume a proportionality to the broadening so that :
\begin{equation}
\begin{aligned}
	T_c(H) - T_{c0}(H) & = \frac{\sigma}{\sigma_0} \left( T_c(0) - T_{c0}(0) \right) \\
	\deltaCT(T_c) & = \deltaCT(T_{c0})\frac{T_c/T_{c0}}{1 + \frac{\sigma_0}{\sigma} \frac{ T_c - T_{c0} }{T_{c0}(0)}}
\end{aligned}
\label{equS:dCT}
\end{equation}
In the last expression, we wrote $T_c = T_c(H)$ and $T_{c0} = T_{c0}(H)$.
As regards the slope, similarly, it should also depend on $T_c$.
Indeed, in high fields for example, where the temperature dependence of $C/T$ is close to linear, the slope should depend both on $T_c$ and on $\deltaCT$.
One way to keep some consistency within the transition is to assume that we have the same entropy balance for all the curves at different $T_c$ at a given field.
At low field, where $C/T(T)$ has no specific reason to remain close to linear far below \Tc, there is no peculiar constraint on this entropy balance (the linear behaviour of $C/T$ below $T_c$ is valid only close enough to $T_c$).
However, for fields closer to $H_{c2}(0)$, we can expect that this entropy balance should be more or less close to zero.
Explicitly, we can enforce that :
\begin{equation}
\begin{aligned}
	\Delta S(T_c) & = \int_0^{T_c} \left[ \deltaCT(T_c) + \alpha(T_c) \, \left( T - T_c \right) \right] dT = \beta T_c \quad \text{, with $\beta$ independent of $T_c$} \\
	\alpha(T_c) &= \frac{2}{T_c} \left[\deltaCT(T_c) - \beta \right] \\
\label{equS:alpha}
\end{aligned}
\end{equation}

Inserting equations (\ref{equS:dCT}) and (\ref{equS:alpha}) in the  equations (\ref{equS:CpIdeal}-\ref{equS:CpDistrib}), we obtain a final expression for $C/T(T)$, easily managed in its integral form by numerical calculations.
It depends linearly on the parameters $\gamma$, $\deltaCT(T_{c0})$ and $\beta$ (close to zero in high fields), and non linearly on $\sigma$ and on $T_{c0}$.
It has two additional inputs, taken from the zero field data: $\sigma_0$ and $T_{c0}(0)$.

\clearpage
\subsection{Measurements of \texorpdfstring{\Hc}{Hc2} for \texorpdfstring{$H \parallel a$\xspace}{H parallel to a}}

Three different samples have been measured and their \Hc determined from the specific heat anomaly for $H \parallel a$\xspace including at very low fields.
The three samples come from different batches. Sample \diese1 and 2 have been measured with the same set up, and sample \diese5 with a different one.
They all exhibit a strong negative curvature near \Tc, proving that this feature is reproducible and intrinsic.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS9_Curvature_H_a.pdf}
	\caption{\Hc of 3 different samples for $H \parallel a$\xspace: the anomalous strong curvature near \Tc is reproducible, even for crystals with very different \Tc.}
	\label{figS:Curvature_H_a}
\end{figure}

\clearpage
\subsection{Measurements of \texorpdfstring{\Hc}{Hc2} for \texorpdfstring{$H \parallel b$\xspace}{H parallel to b} - comparison with resistivity.}

Resistivity has also been measured on the sample from which we cut off sample \diese3.
The critical field obtained with $R=0$ as criterion, is compared to \Hc determined by specific heat.
For the \lftransi phase, as expected, $R=0$ is above the specific heat transition.
At low fields, there is a large difference between the initial slopes at \Tc for the determination from resistivity or specific heat anomaly, most likely due to the sensitivity of resistivity measurements to
filamentary superconductivity, rapidly suppressed by (small) magnetic fields.
For the \hftransi phase, $R=0$ is below the maximum of the specific heat transition, which is more unusual.
This can be seen on Fig.~\ref{figS:RawDataCompCpRho}, showing the temperature dependence of $C/T$ and of the resistivity at a fixed field of 18~T.

This discrepancy can arise from extrinsic inhomogeneities, like a continuous gradient of \Hm in the sample, or from more intrinsic phenomena like a weaker pinning of vortices in the \hftransi phase, which would induce a smaller critical current and possibly a shift of the resistive transition to lower temperatures.
This is well known in organic superconductors \cite{BelinJSuper1999} or in High-Tc cuprates \cite{GrissonnancheNatCom2014}, where the resistivity remains non-zero in the vortex liquid state, favoured by the highly 2D anisotropy of their normal and superconducting properties.
It has also been observed in iron-based superconductors \cite{KoshelevPRB2019}, and like in the organics or high-Tc cuprates, with much stronger differences on the \Tc determination than observed in \UTe.
The difficulty for such an explanation in \UTe is the same as faced for \ucoge \cite{Wu2018}:
the systems are 3D rather than 2D, hence superconducting fluctuations should be much less effective.
In addition, the discrepancy between resistivity and specific heat determination of \Hc occurs only at very high fields, whereas in the other systems, it arises very fast when entering the mixed state.
It could be that the quantitative difference in the effect arises precisely because superconducting fluctuations are much less important in \UTe or \ucoge than in the quasi 2D systems.
However, it remains to be explained why this would happen only in the field-reinforced phase.
This point remains a fully open question.

\begin{figure}[h]
   \begin{minipage}{0.48\textwidth}
     \centering
	\includegraphics[width=1\linewidth]{PicturesSupp/FigS10_PhaseDiagCpRho.pdf}
	\caption{\Hc for $H \parallel b$\xspace determined by $C/T$ measurements (see main article).
	Gray points correspond to the $R=0$ determined by resistivity measurements.
	The shaded region indicates the width of the \hftransi transition.
	The inset is a zoom for field below 4~T.
	}
	\label{figS:PhaseDiagCpRho}
   \end{minipage}\hfill
   \begin{minipage}{0.48\textwidth}
     \centering
	\includegraphics[width=1\linewidth]{PicturesSupp/FigS11_RawDataCompCpRho.pdf}
	\caption{$C/T$ as function of temperature compared to $R$ as function of temperature, both measured at 18T. }
	\label{figS:RawDataCompCpRho}
   \end{minipage}
\end{figure}

\clearpage
\subsection{Comparison of \texorpdfstring{\Hc}{Hc2} and \texorpdfstring{\Hcl}{Hc1}}

Figure~\ref{figS:CompHc2_Hc1} shows \Hc determined by specific heat on sample \diese2 and the \dHc at \Tc calculated and rescaled from \Hcl measurements in ref.~\cite{Paulsen2021}.
The measurements of \Hcl have been done on a sample from the same batch as sample \diese1, and have roughly the same \Tc at zero field.
Near \Tc, in the Ginzburg-Landau regime, the Ginzburg-Landau parameter $\kappa$ is determined by $H_{c1}=\frac{H_c}{\sqrt{2} \kappa} \left(ln(\kappa)+0.49\right)$.
\dHc can then be deduced form  $H_{c2}=\sqrt{2} \kappa H_c$, in order to evaluate \dHc from \dHc
The thermodynamic critical field $H_c(T)$ is determined by double integration of the specific heat at 0~T, and we obtain  $\frac{dH_c}{dT_c}$=-0.056~T/K.
Then, the calculated \dHc are rescaled by the ratio of the respective critical temperatures of sample \diese2 and \diese1 in order to obtain  the values of \dHc for each axis.
These values are displayed in table~\ref{tabS:ParameterTable}.

\begin{figure}[h]
	\centering
	\includegraphics[width=.5\linewidth]{PicturesSupp/FigS12_CompHc2_Hc1.pdf}
	\caption{\Hc measurements on sample \diese2: lines show \dHc calculated and rescaled from \dHcl determined in ref~\cite{Paulsen2021}. The insert is a zoom for fields $H \parallel a$\xspace below 0.25~T. There is a very good agreement between \Hcl and \Hc for $H \parallel a$\xspace and $H \parallel c$\xspace, but a strong discrepancy for $H \parallel b$\xspace.}
	\label{figS:CompHc2_Hc1}
\end{figure}

\begin{table}[h]
    \centering
    \begin{tabular}{c|c|c|c|c|c}
 & -\dHcl (T/K) - ref~\cite{Paulsen2021} & $\kappa$  & -\dHc (T/K) &  -\dHc rescaled (T/K)  &  -\dHc (T/K) - \diese2\\
 \hline
$H \parallel a$\xspace & 0.00113 & 202.683 & 16.052 & 20.480 & 20\\
$H \parallel b$\xspace & 0.00227 & 86.482 & 6.849 & 8.738 & 34.5\\
$H \parallel c$\xspace & 0.00252 & 75.838 & 6.006 & 7.663 & 7.5 \\
\hline
    \end{tabular}
    \caption{Table with the values of the slope \dHcl of \Hcl at \Tc from \cite{Paulsen2021}, the corresponding value of the calculated Ginzburg-Landau parameter $\kappa$, and the predicted value for \dHc for the sample of Ref.~\cite{Paulsen2021} (same batch as \diese1). For sample \diese2, \dHc is rescaled by the ratio of the \Tc of these samples. Last column is the initial slope measured on sample \diese2.}
    \label{tabS:ParameterTable}
\end{table}

\clearpage
\subsection{Field-dependence of the pairing strength modelled by a strong coupling parameter \texorpdfstring{\llambda}{lambda}}

Figure~\ref{figS:CompHc2_BW} shows \Hc along the 3 crystallographic directions, calculated with a strong coupling model for the upper critical field already used in ref.~\cite{Wu2017,Knebel2019}, at fixed pairing strength.
The measured initial slopes \dHc at \Tc are used to determine \vF.
The strong coupling constant $\lambda$ is set to 1, which seems a reasonable value for \UTe.
The plain lines in Fig.~\ref{figS:CompHc2_BW} are \HcT calculated with the orbital limit adjusted to match the measured \dHc, and the gyromagnetic factor $g$ adjusted to match the initial negative curvature along the $a$ and $b$ axes ($g=6.5$ along the \aaxis and $g=0.8$ along the \baxis).
Deviations of the measured \Hc to such a usual combined orbital and paramagnetic limitation are observed for all applied directions of the magnetic field.

\begin{figure}[h]
	\centering
	\includegraphics[width=.5\linewidth]{PicturesSupp/FigS13_CompHc2_BW.pdf}
	\caption{\Hc determined by specific heat on sample \diese2.
	The lines are the best adjustment of \Hc (orbital and paramagnetic limitations) to match the measured initial slopes and curvatures in each direction.
	The dashed doted line corresponds to a pure orbital limitation of \Hc adjusted on its initial slope  for $H \parallel a$\xspace, evidencing the very strong negative curvature close to \Tc.
    $g$ is the gyromagnetic factor and $\lambda$ the strong coupling constant.}
	\label{figS:CompHc2_BW}
\end{figure}

To explain these deviations, a field dependent pairing strength is assumed.
It is extracted from the data through a calculation of \HcT at fixed values of \llambda (see Fig.~\ref{figS:FitsHc2LambdaFixed} for the case of $H \parallel b$\xspace). The typical energy controlling \Tc ($\Omega$), the Coulomb repulsion parameter $\mu^{\ast}$, and the bare average Fermi velocity for field along the $i$ axis ($\bar{v}_{F0}^{bare,i}$) controlling the orbital limit are taken independent of \llambda.
The effective Fermi velocity controlling the orbital limit and so \dHc (at fixed \llambda)  is renormalised as $\bar{v}_{F}^i = \frac{\bar{v}_{F0}^{bare,i}}{1+\lambda}$.
If \llambda is field dependent, this effective Fermi velocity is also field dependent.
In the table~\ref{tabS:ParameterFitHc2}, we also calculate the Fermi velocity along each $i$-axis: $v_{F}^i$, deduced from the effective Fermi velocities through $v_{F}^i = \frac{\bar{v}_{F}^j\bar{v}_{F}^k}{\bar{v}_{F}^i}$, where $j$, $k$ are the axis perpendicular to $i$.

\begin{table}[h]
    \centering
    \begin{tabular}{c|c|c|c|c|c}
 & $\bar{v}_{F}^i(H=0)$ (m/sec) & $g$  & $\bar{\xi}_0^i$ ($\mathring{A}$) &  $v_{F}^i(H=0)$ -along $i$ axis  &  $\xi_0^i$ ($\mathring{A}$)
 - $i$ along axis\\
 \hline
$H \parallel a$\xspace             & 5400 & 0 & 40 & 14400 & 106\\
$H \parallel b$\xspace(\lftransi)  & 8600 & 0 & 64 & 5680 & 42\\
$H \parallel c$\xspace             & 9044 & 0 & 67 & 5130 & 38\\
$H \parallel b$\xspace(\hftransi)  & 8600 & 2 & 64 & 5680 & 42\\
\hline
    \end{tabular}
    \caption{Parameter values of the fit.
    We used a strong-coupling parameter $\lambda(H=0)=1$, with a typical energy (equivalent to the Debye energy) $\Omega=28.4$~K, $\mu^\ast=0.1$, pair-breaking impurity scattering rate $\Gamma=1.39$~K.
    Values of $\bar{v}_{F}$ used in the fit are reported for each field direction. Difference between $\bar{v}_{F}^i$ and $v_{F}^i$ is explained in the text. The corresponding coherence length are calculated from and $\xi_0 = 0.18 \frac{\hbar v_F}{k_B T_c}$ }
    \label{tabS:ParameterFitHc2}
\end{table}

For $H \parallel a$\xspace, as discussed in the main text, we have chosen the most natural hypothesis of an ESP (equal-spin-pairing) state along the \aaxis, hence no paramagnetic limitation of \Hc, an initial slope matching the measured one (agreement with \Hcl) implying a suppression of the pairing strength under field.
However, it is also possible to construct a model where the pairing strength would increase along the \aaxis.
Indeed, maintaining an initial slope matching the measured one, if we suppose that the g-factor is in reality $\ge 6.8$, the same fitting procedure will lead to a field increase of $\lambda(H)$.
We rejected this scenario due to the very large and rather unrealistic value required for the g-factor, and the strong contradiction with NMR Knight-shift measurements observing no change a all for $H \parallel a$\xspace \cite{FujibayashiJPSJ2022}.

\begin{figure}[ht]
	\centering
	\includegraphics[width=.5\linewidth]{PicturesSupp/FigS14_FitsHc2LambdaFixed.pdf}
	\caption{Crosses and squares are data of \Hc for $H \parallel b$\xspace.
	Dashed lines are the \Hc calculated for different fixed values of the coupling constant $\lambda$.
	The green lines correspond to calculation for $g=0$ by steps of $\Delta \lambda = 0.05$, and the blue lines for $g=2$ by steps of $\Delta \lambda = 0.2$.}
	\label{figS:FitsHc2LambdaFixed}
\end{figure}

\clearpage
\subsection{Broadening of the specific heat transition by a distribution of \texorpdfstring{\Hm}{Hm}}
As explained in section \ref{secS:SuppCpHm} of these supplement, we could extract a standard deviation $\sigma = 0.19 $~T for the distribution of \Hm (see Fig.~\ref{figS:param_transi_Hm}), hence a relative standard deviation $\frac{\sigma}{H_m} \sim 0.55 \%$

As explained in the main text, from the calculation of \Hc at fixed values of the pairing strength \llambda (Fig.\ref{figS:FitsHc2LambdaFixed}), we can extract also the superconducting critical temperature under field as:
\begin{equation}
    \begin{aligned}
        T_c & = \varphi\left( H, \Tilde{\lambda} \left( \frac{H}{H_m} \right) \right) \\
        \Tilde{\lambda} \left( \frac{H}{H_m} \right) & = \lambda \left( H \frac{H_{m0}}{H_m} \right)
    \end{aligned}
    \label{equS:Tc_H_lambda}
\end{equation}
Where $H_{m0}$ is the centre of the distribution of metamagnetic fields \Hm, determined from the specific heat measurements of the metamagnetic transition.
\llambdaH is the field dependent pairing strength deduced from the different models for \Hc and drawn on Fig.~\ref{fig:Lambda_H} of the main text.

From this relation, we can calculate the effect of a Gaussian distribution of \Hm on the specific heat anomaly at constant field of the superconducting transition, using for \CT (instead of Equ.~\ref{equS:CpDistrib}):
\begin{equation}
	C/T = \int p(H_m)  C/T\left(T,T_c\left(H,H_m\right)\right)  dH_m
	\label{equS:Cp_distrib_Hm}
\end{equation}
This is the way we could draw the broadening of the specific heat anomaly on Fig.~\ref{fig:FitCpHF} and Fig.~\ref{fig:FitCp_theta} of the main text, using the two different determinations of \llambdaHHm (with or without paramagnetic limitation of \Hc).

However, even without a full determination of the shape of the anomaly, requiring a numerical integration of Equ.~\ref{equS:Cp_distrib_Hm}, we can understand why the broadening is larger when there is a paramagnetic limitation of \Hc.
From Equ.~\ref{equS:Tc_H_lambda}, we can derive the derivative of \Tc with respect to \Hm at fixed H and for $H_m=H_{m0}$.
It measures the sensitivity of \Tc to \Hm, hence the broadening of the \CT anomaly due to a distribution of \Hm:
\begin{equation}
    \left.\frac{\partial T_c}{\partial H_m}\right\vert_H =  \left.\frac{\partial T_c}{\partial \lambda}\right\vert_H \left(-\frac{H}{H_{m0}}\right) \left(\frac{d \lambda}{d H} \right)
    \label{equS:dTc_Hm}
\end{equation}
When comparing both models, it is clear that one has a stronger field dependence of \llambda than the other, but this could be compensated by a different $\left.\frac{\partial T_c}{\partial \lambda}\right\vert_H$ which has to be computed at finite field (on the \HcT line).
Indeed, both models share the same \HcT.
We can compute its temperature derivative from Equ.~\ref{equS:Tc_H_lambda}:
\begin{equation}
    \begin{aligned}
        d T & =  \left.\frac{\partial T_c}{\partial H}\right\vert_{\lambda} dH_{c2} +
        \left.\frac{\partial T_c}{\partial \lambda}\right\vert_{H} \left(\frac{H_{m0}}{H_{m}}\right) \left(\frac{d \lambda}{d H} \right)dH_{c2} \\
        \frac{dT}{d H_{c2}} - \left.\frac{\partial T_c}{\partial H}\right\vert_{\lambda} & =  -  \left(\frac{H_{m0}}{H}\right) \left.\frac{\partial T_c}{\partial H_m}\right\vert_H
    \end{aligned}
    \label{equS:dTc_Hc2}
\end{equation}
The last equation shows that the difference between models for $\left.\frac{\partial T_c}{\partial H_m}\right\vert_H $ arises not directly from $\left(\frac{d \lambda}{d H} \right)$, but rather from  $\left.\frac{\partial T_c}{\partial H}\right\vert_{\lambda}$:
this term is much larger when \HcT at fixed \llambda becomes "flat" due to the paramagnetic limitation arising for singlet pairing, than for a pure orbital limit in case of a spin-triplet ESP state (see Fig.~\ref{figS:FitsHc2LambdaFixed}).

\clearpage
\subsection{Angular dependence in the \texorpdfstring{($b$,$c$)}{(b,c)} plane}
AC specific heat measurements have been performed on sample \diese2 up to 18.5 T for several angles in the ($b$,$c$) plane.
Figure~\ref{figS:Cp_vs_H_3_angles} shows temperature sweeps for different fields for angles of 0$^\circ$, 10$^\circ$ and 15$^\circ$ from $b$ toward the $c$ axis.
As the field is rotated toward the $c$ axis, the sharp transition of the \lftransi phase is shifted toward lower temperatures.
The same behaviour is observed for the wide transition of the \hftransi phase.
The corresponding critical temperature for the two transitions at the different angles are reported on the phase diagram of Fig.~\ref{figS:PhaseDiag_vs_angle}, using the same Gaussian analysis as in the main paper.
For \Hc at 15$^\circ$, it was impossible to extract a reliable value of \Tc for the \hftransi transition.

\begin{figure}[b]
     \centering
	\includegraphics[width=0.31\textwidth]{PicturesSupp/FigS15a_Cp_vs_H_0degree.pdf}
		\hspace{\fill}
	\includegraphics[width=0.31\textwidth]{PicturesSupp/FigS15b_Cp_vs_H_10degree.pdf}
		\hspace{\fill}
	\includegraphics[width=0.31\textwidth]{PicturesSupp/FigS15c_Cp_vs_H_15degree.pdf}
\caption{$C/T$ as a function of temperature for several fields between 0~T and 18.5~T, measured on sample \diese3.
(a) $C/T$ for $H \parallel b$\xspace,
(b) for an angle of 10$^\circ$ toward the \caxis and (c) for an angle of 15$^\circ$ toward the \caxis.}
\label{figS:Cp_vs_H_3_angles}
\end{figure}

\begin{figure}[hb]
	\centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS16_PhaseDiag_vs_angle.pdf}
	\caption{Phase diagram for a field applied along the \baxis (red squares and crosses), at 10$^\circ$ toward $c$ axis (black triangles and circles), and at 15$^\circ$ (orange diamonds). The empty marker represents fields where the width and jump sizes of the \hftransi transitions need to be fixed to determine \Tc from a Gaussian fit of both anomalies.}
	\label{figS:PhaseDiag_vs_angle}
\end{figure}

Similarly, as presented in the main article in Fig.~\ref{fig:FitCp_theta} for an applied field of 18.5~T, temperature sweeps have been performed at 16~T for different angles in the ($b$,$c$) plane. Results are presented in Fig.~\ref{figS:Cp_16T_vs_angle}. The \lftransi and \hftransi transitions overlap and the \hftransi is less defined as at 18.5~T. Nevertheless, we clearly see a shift to the low temperatures of the \lftransi transition as the field is rotated away from \baxis. The \hftransi transition also clearly disappears, shifting to the low temperatures and the jump size decreasing.

\begin{figure}
	\centering
	\includegraphics[scale=0.5]{PicturesSupp/FigS17_Cp_16T_vs_angle.pdf}
	\caption{$C/T$ as a function of temperature at 16 T, for different angles in the ($b$,$c$) plane.}
	\label{figS:Cp_16T_vs_angle}
\end{figure}

\clearpage