\section*{Code Metadata}
    \noindent Current code version: v0.5.0 \\
    Link to repository: \url{https://github.com/CenterForTheBuiltEnvironment/clima} \\
    Legal Code License: MIT \\
    Code versioning system used: git \\
    Software code language: Python \\
    Documentation: \url{https://center-for-the-built-environment.gitbook.io/clima/} \\
    Support email: \url{https://github.com/CenterForTheBuiltEnvironment/clima/issues}\\

    \section*{Nomenclature}

\renewcommand{\baselinestretch}{0.75}\normalsize
\renewcommand{\aclabelfont}[1]{\textsc{\acsfont{#1}}}
\begin{acronym}[longest]

    \acro{t-db}[$T_{db}$]{dry-bulb air temperature\acroextra{, $^{\circ}$C}}
    \acro{t-wb}[$T_{wb}$]{wet-bulb air temperature\acroextra{, $^{\circ}$C}}
    \acro{t-dp}[$T_{dp}$]{dew point air temperature\acroextra{, $^{\circ}$C}}
    \acro{ghi}[$GHI$]{Global Horizontal Irradiance\acroextra{, Wh/m$^2$}}
    \acro{ws}[$U$]{wind speed\acroextra{, m/s}}
    \acro{hr}[$HR$]{Humidity Ratio\acroextra{, kg$_{water}$/kg$_{dry\;air}$}}
    \acro{t-op}[$T_{o}$]{operative air temperature\acroextra{, $^{\circ}$C}}
    \acro{tg}[$t_{g}$]{globe temperature\acroextra{, $^{\circ}$C}}
    \acro{rh}[$RH$]{relative humidity\acroextra{, \%}}
    \acro{v}[$V$]{average air speed\acroextra{, m/s}}
    \acro{t-r}[$\overline{t_{r}}$]{mean radiant temperature\acroextra{, $^{\circ}$C}}
    \acro{clo}[$I_{cl}$]{total clothing insulation\acroextra{, clo}}
    \acro{i-cl}[$i_{cl}$]{permeation efficiency of water vapor through the clothing layer}
    \acro{met}[$M$]{rate of metabolic heat production\acroextra{, W/m\textsuperscript{2}}}

    \acro{e}[$\varepsilon$]{average emissivity of clothing or body surface}
    \acro{sigma}[$\sigma$]{Stefan-Boltzmann constant\acroextra{, 5.67 x 10\textsuperscript{-8} W/(m\textsuperscript{2}K\textsuperscript{2})}}

    \acro{pmv}[PMV]{Predicted Mean Vote}
    \acro{ppd}[PPD]{Predicted Percentage of Dissatisfied\acroextra{, \%}}
    \acro{set}[SET]{Standard Effective Temperature\acroextra{, $^{\circ}$C}}
    \acro{ce}[CE]{Cooling Effect\acroextra{, $^{\circ}$C}}
    \acro{phs}[PHS]{Predicted Heat Strain}

    \acro{BMS}[BMS]{Building Management System}
    \acro{HVAC}[HVAC]{Heating, Ventilation, and Air Conditioning}
    \acro{VAV}[VAV]{Variable Air Volume}
    \acro{AHU}[AHU]{Air Handling Unit}

    \acro{ema}[EMA]{Ecological Momentary Assessment}
    \acro{sdk}[SDK]{Software Development Kit}
    \acro{api}[API]{Application Programming Interface}

    \acro{wmo}[WMO]{World Meteorological Organization}
    \acro{who}[WHO]{World Health Organization}
    \acro{cdc}[CDC]{Centers for Disease Control and Prevention}
    \acro{noaa}[NOAA]{National Oceanic and Atmospheric Administration}
    \acro{epa}[EPA]{United States Environmental Protection Agency}
    \acro{iea}[IEA]{International Energy Agency}
    \acro{un}[UN]{United Nations}
    \acro{epw}[EPW]{EnergyPlus Weather}
    \acro{utci}[UTCI]{Universal Thermal Climate Index\acroextra{, $^{\circ}$C}}
    \acro{svg}[SVG]{Scalable Vector Graphics}
    \acro{csv}[CSV]{Comma-Separated Values}

\end{acronym}
\renewcommand{\baselinestretch}{1}\normalsize

    \section{Motivation and significance}\label{sec:motivation-and-significance}
Buildings construction and operation account for $\approx$~36\% of global final energy use and nearly 40\% of energy-related carbon dioxide (CO\textsubscript{2}) emissions~\cite{InternationalEnergyAgency2018}.
However, a sixth of the occupants of the building are not satisfied with their indoor environment~\cite{Graham2021}.
There is a need to reduce the greenhouse gas emissions related to buildings and increase their resilience while improving occupant well-being.
Energy consumption can be reduced and occupant comfort can be improved when a building is climatically adapted to its environment~\cite{day_gunderson_2015, brager_baker_2009, altomonte_schiavon_kent_brager_2017}.
The climatically adapted design of a building must take into account the outdoor and indoor climates, optimize site selection, building orientation relative to the sun and the wind, daylight, and view access~\cite{SANTAMOURIS201374, RATTI200349}.
This process is present in vernacular design traditions all over the world~\cite{mileto2014vernacular}.
The early modern movement in architecture took a renewed interest in questions related to solar exposure and orientation, as can be clearly recognized in many sketches, built works, and theoretical writings by Le Corbusier and Alvar Aalto~\cite{menin2003nature}.
After World War II, data and physics-based approaches to buildings climate adaptation started to emerge.
In the early 1950s, the Form and Climate Research Group at the Columbia Graduate School of Architecture sought to develop techniques for refining the design process according to climatic adaptability~\cite{Barber2020MAaC}.
Almost contemporary, the American Institute of Architects started publishing a seminal set of climate graphs to describe the climate of different American cities~\cite{Taylor1949}.
This work about climate data, climatic conditions, and their implication for architectural design was then brought forward by the research of Victor and Aladar Olgyay, among whose many contributions there is the Bioclimatic Chart~\cite{Olgyay1963}.
Shortly after, the emergence of Computer-Aided Design generated a few early tools that aimed at creating climatic data visualizations for architects and building engineers, such as the work by Murray Milne at UCLA~\cite{Milne}.
From 1978 onward, Milne has developed several tools that use climatic data as design guidance, most notably Climate Consultant~\cite{climateconsultant}.
One of the early challenges for the Architecture Engineering and Construction (AEC) industry was the lack of a consistent climate data format.
This issue was solved with the introduction of the \ac{epw} format~\cite{Drury}.
In the early 2000s, Andrew Marsh developed Ecotect and the Weather Tool, now retired, which have enabled thousands of practitioners and students to integrate climate data in their CAD workflows for the first time~\cite{ECOTECT}.
The Weather Tool shared with Climate Consultant the fact of being a stand-alone, easy-to-use application that enabled the visualization of weather data with some degree of customization.
Recently, new tools have emerged in the field of climate data visualization for the AEC industry, most notably Ladybug Tools~\cite{Ladybug}, Climate Studio~\cite{ClimateS53:online}, and cove.tool~\cite{covetool38:online}.
All of these projects enable different climate visualization directly in the most prevalent CAD environments, mainly Autodesk Revit and Rhinoceros 3D~\cite{mcneel-rhinoceros}.
To achieve this integration, the tools act as software plug-ins.
This strategy has certain advantages, i.e., showing contextual information in the user’s design software of choice.
On the other hand, it makes the tools part of a closed-source ecosystem and creates barriers to their use as they require multiple installations and increasing levels of specialization.
For example, to use Ladybug one needs to be familiar with Grasshopper, the visual programming extension to Rhinoceros 3D~\cite{mcneel-rhinoceros}.
Existing stand-alone tools, like Climate Consultant, do not provide vector graphics and therefore their use in practice and education is limited~\cite{climateanalysistool}.

Our aim is to develop open-source software that overcomes the limitations listed above, reduce barriers to access, and helps users to design and construct healthier and sustainable buildings.

\section{Software description}\label{sec:software-description}
The \gls{clima} is a free open-source web application for the analysis and visualization of climate data specifically designed to support the needs of architects, engineers, students, and educators interested in climate-adapted building design.
It does not only visualize data, but it organizes, manipulates, and relates them in an easy-to-interpret fashion.
The landing page is characterized by an interactive world map.
Through this map, the user can access circa 30,000 publicly available weather data files.
Data are collected from publicly available repositories maintained by EnergyPlus and Climate.OneBuilding.org.
Alternatively, the user can upload any valid \ac{epw} file.
The uploaded file is automatically processed and analyzed.
Our tool does not perform quality checks other than checking that the data supplied matches the expected structure of an \ac{epw} file.
In case of a mismatch, an error message is shown to the user.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\columnwidth]{figures/thumbnail.png}
    \caption{Video tutorial -- An introduction to the CBE Clima Tool.}
    \label{fig:video}
\end{figure}

The data contained within the \ac{epw} files are used to generate easy-to-interpret charts.
Many of these allow the user to explore data at different temporal resolutions, highlighting seasonal, daily, or hourly patterns.
All the charts allow some degree of interactivity, allowing the user to hover, zoom, and highlight data.
For some charts, we also allow the user to select specific time ranges (i.e. wind roses analysis).
Not all relevant data for a climate analysis is contained directly in the \ac{epw} files but can be derived from the information contained therein.
We use the open-source Python library \verb|pvlib|~\cite{F.Holmgren2018} to calculate the apparent solar position (altitude and azimuth).
We use the open-source Python library \verb|pythermalcomfort|~\cite{pythermalcomfort} to calculate other indices that can be derived from those contained in the EPW files.
For example, we calculate the \ac{utci}~\cite{utci} scenarios to understand outdoor comfort, solar gain on occupants~\cite{Arens2015}, the adaptive thermal comfort model~\cite{ASHRAE552020}, natural ventilation potential, heating and cooling degree days, and various other psychrometric variables.
The \gls{clima} can be accessed at \url{https://clima.cbe.berkeley.edu}.

\subsection{Software Architecture}\label{subsec:software-architecture}
The \gls{clima} tool is built using Dash~\cite{plotly} a Python framework to build web analytic applications and data visualization websites.
We use Dash since it is an open-source library released under an MIT license, and applications built with it can be rendered in any web browser.
Clima is, therefore, cross-platform and compatible with the most widely used web browsers.
Dash applications are web servers that run Flask and communicate JSON packets over HTTP requests.
We have deployed \gls{clima} using a cloud service that hosts our containerized application.
We use Docker to create a container that listens to user requests.
A request is triggered each time the user navigates to the URL: \url{https://clima.cbe.berkeley.edu}.
Dash’s front-end renders components using React.js, a Javascript framework.
While the back-end code runs on a cloud instance, the front-end code is run locally inside the user's browser.
Each time a user changes some of the input parameters, a new request is sent to the server.
The Python source code is stored inside the container, while the EPW data is queried from the databases maintained by EnergyPlus and Climate.OneBuilding.org.
A diagram depicting the application architecture is shown in Figure~\ref{fig:application_arch}.
For more information about how a Dash application works please visit and consult their official documentation.

\begin{figure*}[htb!]
    \centering
    \includegraphics[width=\textwidth]{figures/applicationArchitecture}
    \caption{Application architecture diagram.}
    \label{fig:application_arch}
\end{figure*}

\subsection{Software Functionalities}\label{subsec:software-functionalities}
The \gls{clima} allows users to analyze and visualize climatic data.
Users can access \ac{epw} files via a map-based interface (as previously stated in Section~\ref{sec:software-description}).
All images and charts included in our tool can be downloaded in \ac{svg} format.
This allows them to be modified and adapted according to the users' needs.
Furthermore, the vector nature of the \ac{svg} format enables both high-resolution on-screen display and high-quality printing.
This is an important feature of the tool as it enables visually compelling data communication with multiple stakeholders.
The interactive nature of the graphs means that the user can inspect the underlying numerical information through rich mouse-over events.
Users can download, in \ac{csv} format, either the \ac{epw} file or the \gls{clima} data frame created in the app back-end to generate the plots.
The latter contains the variables in the \ac{epw} and additional ones derived from it, as previously described.
When displaying the variables in the interactive plots, users can either select to use the Global or Local Value Ranges, using the radio button at the top of the page, see Figure~\ref{fig:clima_home}.
The `Global' option uses preset limits for the chart axes, chosen to cover the vast majority of the climatic ranges to be found on Earth.
This allows users to effectively compare charts generated for different locations.
The `Local' option sets the upper and lower limits of the chart axes as a function of the data contained in the \ac{epw} file.

The following paragraphs give an overview of the main functionalities of each page of the tool.
However, we recommend consulting the official project documentation for more detailed and up-to-date information available at \url{https://cbe-berkeley.gitbook.io/clima/}.

\paragraph{Home Page} Users can either choose to analyze the climate of the locations displayed on the map or upload a custom \ac{epw} file as shown in Figure~\ref{fig:clima_home}.
After loading an \ac{epw} file the user can then access the other tabs.

\begin{figure*}[htb!]
    \centering
    \includegraphics[width=\textwidth]{figures/clima_home}
    \caption{\gls{clima} home page where the \ac{epw} file is selected.}
    \label{fig:clima_home}
\end{figure*}

\paragraph{Climate Summary} It provides textual information about the selected location such as longitude, latitude, elevation, Koppen-Geiger climate zone, average/hottest/coldest yearly temperatures, annual cumulative horizontal solar radiation, and percentage of diffuse horizontal solar radiation.
The user is presented with a further map showing the location of the weather station selected.
The user can also download the source data, visualize the heating and cooling degree day chart and compute them as a function of user-defined base temperatures and analyze the distributions of the \ac{t-db}, \ac{rh}, \ac{ghi}, and \ac{ws}.

\paragraph{Temperature and Humidity} It shows the yearly, hourly, and heat-map charts for \ac{t-db} and \ac{rh}.
The yearly \ac{t-db} chart depicts the daily temperature range and the mean value together with the ASHRAE adaptive thermal comfort region for the 80\% and 90\% acceptability ranges.
The prevailing mean outdoor temperature was calculated using an exponentially weighted running mean of the daily outdoor temperature for the last 7 days and assuming an $\alpha$ of 0.9~\cite{EN1679812019}.
The daily chart shows how the hourly readings vary as a function of the hour of the day, and data are grouped by month.
The heat map shows each individual hourly data point in a grid where the vertical axis shows the time of the day, while the horizontal axis shows the day of the year.
The table with descriptive statistics for each month and year is shown at the bottom of the page.

\paragraph{Sun and Clouds} It visualizes the apparent sun path as seen by an observer at the weather station location either in Spherical coordinates (where the observer sits at the weather station location and the apparent altitude and azimuth are mapped to the polar and azimuthal angles respectively) or Cartesian coordinates (where the apparent altitude and azimuth are mapped to a Cartesian plot).
The monthly average hourly global and diffuse solar radiation (kWh/m\textsuperscript{2}), and the monthly cloud coverage are also presented.
Users can also plot and visualize how any of the many variables related to solar radiation contained in the \ac{epw} file (e.g. \ac{ghi}, extraterrestrial horizontal irradiation) vary throughout the year.

\paragraph{Wind} It provides information about the wind velocity and direction for the selected location.
It comprises several wind roses that either display the full dataset (yearly data), values grouped by seasons or grouped by time of day.
The user can also create a custom wind rose by filtering the data by month and hour of the day.
This page also shows two heat maps, one for the wind speed and one for the wind direction.

\paragraph{Psychrometric Chart} It shows data on a customizable psychrometric chart.
Users can select to bin the data and colour them by frequency or display the data and colour them by any variable of choice.
Optionally, filters can be applied to select data contained in a specific time range, or filtered according to a secondary variable.
As an example, see the figure in Section~\ref{subsec:psy} for more.

\paragraph{Natural Ventilation} It allows the user to estimate when and for how long the climate is suitable for the use of natural ventilation.
This can be achieved by selecting an outdoor \ac{t-db} range in which natural ventilation can be used.
The user can also filter data either based on monthly or hourly ranges.
If radiant systems are used to cool the building, the user can also specify the surface temperature at which they operate.
The \gls{clima} filters out automatically data when the surface of the radiant system is below the dew point.

\paragraph{Outdoor comfort} It plots the \ac{utci} values on two heat-maps.
The first chart shows the \ac{utci} equivalent temperature, while the second shows the data binned into the 10 \ac{utci} stress categories (from extreme cold stress to extreme hot stress).
Four distinct \ac{utci} scenarios are calculated showing the impact of exposure to the sun and the wind and combinations thereof.
This might help the designer to understand what combination of exposure to or protection from the sun and wind would provide the most comfortable outdoor environment at different times of the year.

\paragraph{Data Explorer} It allows the user to explore and plot all variables contained in the source dataset.
It comprises three charts (yearly chart, daily chart, and heat-map chart) that are also displayed on the `Temperature and Humidity' page.
Moreover, it also allows users to visualize the correlation between three input variables.
For example in Figure~\ref{fig:clima_explore} is a screenshot of the aforementioned charts, and we are showing how data points are distributed as a function of \ac{t-db}, \ac{rh}, and \ac{ghi}.

\begin{figure*}[htb!]
    \centering
    \includegraphics[width=\textwidth]{figures/clima_explore.png}
    \caption{Data explorer tab -- The user can select two variables one to be displayed on the x-axis and one to be displayed on the y-axis.
    This allows the user to visualize the correlation between the two variables.
    While the histograms depict the distribution of the data.
    This selection is then used to create a scatter plot and a heat-map.
    In the scatter plot a third variable is used to color the markers and the range is displayed in the color bar.
    The heat-map, on the other hand, depicts the correlation between the x and y variables and shows the number of points that are comprised in each bin.
    }
    \label{fig:clima_explore}
\end{figure*}

\section{Illustrative examples}\label{sec:illustrative-examples}

\subsection{Sun Path Chart}
Sun control and shading systems can significantly reduce peak heat gain and cooling requirements, while also improving the quality of natural lighting in the interior of the building.
This example shows how the user can visualize how \ac{t-db} varies as a function of the position of the sun in the sky.
This information can then be used to optimize the design of sun-shading devices.
This chart can be generated from the `Sun and Clouds` tab.
The user first has to select, using the drop-down, between a spherical or Cartesian coordinate system.
In this example, we use the latter and decided to select \ac{t-db} as a secondary variable to plot from the second dropdown.
Figure~\ref{fig:clima_sun_path} shows that this location is characterized by cold winters and hot summers.
The user can then determine the optimal size and shape of the shading that can be used to minimize solar gains in summer while maximizing passive solar heating in winter.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\columnwidth]{figures/clima_sun_path.png}
    \caption{Cartesian sun path chart.
    It can be used to determine the solar azimuth and altitude for a given location.
    In addition, we allow users to colour the points using a third variable.
    This allows for example to determine how the \ac{t-db} varies as a function of the solar azimuth and altitude throughout the year.}
    \label{fig:clima_sun_path}
\end{figure}

\subsection{Psychrometric Chart}\label{subsec:psy}
This example shows how the user can visualize the distribution of the climate data for a specific location on the psychrometric chart.
We also decided to colour the data points as a function of the \ac{utci} calculated for the condition where the reference person is neither exposed to direct solar radiation nor to wind.
This information can be used to design a semi-outdoor space or to determine when natural ventilation may be used.
This chart can be generated from the `Psychrometric` tab.
Through the drop-down located on the top left side of the page, the user can select the variable `UTCI: no Sun \& no Wind`.
The \gls{clima} automatically generates the chart shown in Figure~\ref{fig:clima_psy}.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\columnwidth]{figures/clima_psy.png}
    \caption{Psychrometric chart.
    It can be used to visualize the relationship between dry air and moisture.
    In addition, we allow users to colour the points using a third variable, in this example, we are showing how the \ac{utci} -- No Sun and Wind -- varies as a function of the dry-bulb air temperature and the humidity ratio.}
    \label{fig:clima_psy}
\end{figure}

\section{Impact}\label{sec:impact}
The \gls{clima} is intended to be used by architects, engineers, researchers, educators, and students.
It can even be used by non-technical users with no analytical or programming skills.
Compared with similar competing solutions for climate analysis, \gls{clima} offers a very complete and flexible feature set that makes it useful for a varied group of users.
The \gls{clima} can be used by students exposed to this topic for the first time and by designers with many years of experience.
The easy output of high-quality vector graphics licensed under Creative Commons Attribution 4.0 International (CC BY 4.0) facilitates its integration into presentation documents.
While alternative solutions require either subscriptions or one or more installs, \gls{clima} can be used via any device equipped with a browser.
Some of the most unique features of the tool, in addition to the completeness, interactivity and customizability of the data visualizations produced, pertain to the availability of complex metrics such as Natural Ventilation potential and the Universal Thermal Climate Index (UTCI) under four distinct scenarios.
We have also open-sourced the code so users around the world can contribute to our project or use our functions in their applications.
The source code is licensed under the MIT license.
This permissive license enables users to share, copy, and redistribute the material in any medium or format and to adapt, remix, transform, and build upon the material for any purpose, even commercially.
These choices were dictated in the hope of promoting the widest possible adoption.
We released the first version of this tool in August 2021.
Since its release, the tool has consistently gathered over 2000 organic unique monthly users from more than 70 nations worldwide.

\section{Conclusions}\label{sec:conclusions}
The \gls{clima} is a free and open-source web application for the analysis and visualization of climate data specifically designed to support the needs of architects, engineers, students, and educators interested in the design of climate-adapted buildings.
It aims to reduce barriers to the implementation of sustainable design.
\gls{clima} processes and visualizes \ac{epw} climate files.
It allows users to perform complex psychrometric data calculations and interactively visualize the results without the need of writing any code or installing specialized software.
It is written in Python.
All generated charts can be downloaded in \ac{svg} format; this allows them to be modified and adapted according to the users' needs.
It is composed of the following main pages: home, climate summary, temperature and humidity, sun and clouds, psychrometric chart, natural ventilation, outdoor comfort, and data explorer. The pages contain interactive graphs and tables.
The tool allows an easy and accessible interpretation of climate data.

\section{Conflict of Interest}\label{sec:conflict-of-interest}
We confirm that there are no known conflicts of interest associated with this publication and that there has not been financial support for this work that could have influenced its outcome.

\section*{Acknowledgements}
We would like to acknowledge the work of the authors who contributed to the development of the \gls{clima}\footnote{\url{https://github.com/CenterForTheBuiltEnvironment/clima/graphs/contributors}}.
This research has been supported by the Center for the Built Environment at the University of California Berkeley and the Republic of Singapore’s National Research Foundation through a grant to the Berkeley Education Alliance for Research in Singapore (BEARS) for the Singapore-Berkeley Building Efficiency and Sustainability in the Tropics (SinBerBEST) Program.

\bibliographystyle{elsarticle-num}


\section{Introduction}

Topological quantum computation based on Majorana bound states is a leading candidate for processing quantum information~\cite{nayak_non-abelian_2008}. The vortex cores of topological superconductors, such as gapped $p$-wave superconductors, host such self-conjugate Majorana bound states at zero energy~\cite{read_paired_2000, ivanov_non-abelian_2001, alicea_new_2012}. The search for these unconventional superconductors has greatly intensified in the past few years as new candidate $p$-wave pairing states are proposed both intrinsically in bulk superconductors~\cite{maeno_superconductivity_1994, maeno_intriguing_2001, mackenzie_superconductivity_2003, mackenzie_even_2017, ran_nearly_2019, li_observation_2019, xu_spin-triplet_2020, aoki_unconventional_2022} and on the surfaces of more conventional superconductors~\cite{fu_superconducting_2008, wang_topological_2015, wu_topological_2016, xu_topological_2016, zhang_observation_2018, machida_zero-energy_2019}.

The Little-Parks effect~\cite{little_observation_1962} originates from the macroscopic quantum coherence of Cooper pairs; due to the quantization of the fluxoid, the resistance of a thin superconducting ring oscillates as a function of the applied magnetic flux. For a conventional $s$-wave superconductor, the periodicity of these Little-Parks oscillations is given by the flux quantum $\Phi_0 = h/2e$, and the minima of the resistance correspond to integer multiples of $\Phi_0$. It has been recognized, however, that unconventional superconductors may exhibit different kinds of Little-Parks oscillations. For example, in gapless superconductors with $d$-wave pairing, the Little-Parks oscillations acquire an enlarged periodicity 2$\Phi_0$~\cite{loder_magnetic_2008,juricic_restoration_2008,zhu_magnetic_2010}, while polycrystalline $p$-wave superconductors have shifted Little-Parks oscillations with the minima of the resistance corresponding to half-integer multiples of $\Phi_0$~\cite{geshkenbein_vortices_1987,xu_spin-triplet_2020,li_observation_2019}. Fractional Little-Parks oscillations with reduced periodicities $\Phi_0 / n$ have also been recently reported both experimentally~\cite{ge_discovery_2022} and theoretically~\cite{halasz_fractional_2021, zhang_higgs-leggett_2022, pan_frustrated_2022}.

For intrinsic spin-triplet $p$-wave superconductors, Majorana bound states have been predicted to emerge in the cores of half-quantum vortices (HQVs) at which half-integer flux quanta $\Phi_0/2$ pierce through the superconductor~\cite{read_paired_2000, ivanov_non-abelian_2001}. In the presence of such HQVs, the Little-Parks oscillations are then expected to possess a distinctive two-peak structure with minima of the resistance at both integer and half-integer multiples of $\Phi_0$~\cite{vakaryuk_effect_2011, yasui_little-parks_2017, cai_magnetoresistance_2022}. At each half-integer minimum, the fluxoid of the superconducting ring is quantized to a half-integer multiple of $\Phi_0$ (meaning that an HQV is bound to the central hole of the ring), while the two peaks around such a minimum correspond to transitions between integer and half-integer fluxoid quantizations. Still, even though this two-peak structure in the Little-Parks oscillations may prove crucial for identifying spin-triplet superconductors, the precise conditions required for its observation are yet to be firmly established. More generally, a rigorous theoretical understanding of the spin-triplet Little-Parks effect could reveal additional signatures of spin-triplet superconductivity and hence provide alternative avenues for detecting this exotic superconducting state on the device scale.

In this work, we employ the Ginzburg-Landau approach to theoretically study the Little-Parks effect in spin-triplet superconducting rings supporting HQVs. We consider both the ``conventional'' Little-Parks oscillations of the critical temperature~\cite{little_observation_1962} and the analogous magnetoresistance oscillations below the critical temperature that result from thermal vortex tunneling~\cite{cai_magnetoresistance_2022, sochnikov_large_2010, sochnikov_oscillatory_2010, cai_unconventional_2013, mills_vortex_2015}. In computing the residual resistance of the superconducting ring below the critical temperature, we do not only focus on the fluxoid ground state~\cite{cai_magnetoresistance_2022, aoyama_little-parks_2022} but also account for the thermally occupied excited states and the thermally activated transitions between them.

We first demonstrate that HQVs are stabilized by an appropriate higher-order term in the Ginzburg-Landau free energy which penalizes the charge supercurrent but not the spin supercurrent. By analyzing the effect of this term on the free energies of the various fluxoid states below the critical temperature, we confirm the presence of a two-peak structure in the magnetoresistance oscillations~\cite{vakaryuk_effect_2011, yasui_little-parks_2017, cai_magnetoresistance_2022} and understand how the separation between the two peaks depends on the temperature. Next, we compute the magnetoresistance oscillations themselves and explicitly quantify the prominence of the characteristic two-peak structure. By identifying a fundamental trade-off between minimizing the width of each peak and maximizing the overall magnitude of the resistance, we provide detailed experimental guidelines for probing the two-peak structure in real candidate materials.

Turning our attention to the ``conventional'' Little-Parks oscillations, we observe two separate critical temperatures with the higher one marking the onset of superconductivity and the lower one separating a single-spin superconducting state above and a spin-triplet superconducting state below. While the two-peak structure is entirely absent from the Little-Parks oscillations of the upper critical temperature, it translates into a two-valley structure for the lower critical temperature, thus providing a further signature of spin-triplet superconductors supporting HQVs.

\section{Ginzburg-Landau theory for a spin-triplet superconductor}

We consider a spin-triplet superconductor with $p_x+ip_y$ pairing symmetry in which spin-orbit coupling energetically favors $(\uparrow\uparrow)$ and $(\downarrow\downarrow)$ Cooper pairs over $(\uparrow\downarrow) + (\downarrow\uparrow)$ Cooper pairs. Such a superconductor can support HQVs around which the superconducting phase of only one type of Cooper pair [either $(\uparrow\uparrow)$ or $(\downarrow\downarrow)$] winds by $2\pi$. Under a magnetic field parallel to the spin quantization axis, the Ginzburg-Landau free energy of such a superconductor is given by~\cite{machida_incomplete_1978}
\begin{widetext}
\begin{eqnarray}\label{eq:FreeEnergy}
F = \frac{T_c F_0} {V_0} \int d^3r\Bigg\{ \sum_{\sigma = \uparrow,\downarrow}\biggl[-\Big\{(1-t)&+&z_\sigma |\grad \times \vec{a}|\Big\}|\psi_\sigma|^2+\frac{|\psi_\sigma|^4}{2}+\xi_0^2\big|(\grad-i\vec{a})\psi_\sigma\big|^2\biggl] \nonumber \\
&+& c|\psi_\uparrow|^2|\psi_\downarrow|^2+\xi_0^2 d \big|(\grad-2i\vec{a})[\psi_\uparrow\psi_\downarrow]\big|^2 +\frac{V_0} {2\mu_0 T_c F_0}(\grad \times \vec{a})^2\Bigg\},
\end{eqnarray}
\end{widetext}
where $t = T/T_c$ is the dimensionless temperature (with $T_c$ being the critical temperature), $\xi_0$ is the zero-temperature coherence length, $\vec{a} = 2\pi \vec{A}/\Phi_0$ is the magnetic vector potential, while $\psi_\uparrow$ and $\psi_\downarrow$ are the superconducting order parameters corresponding to the $(\uparrow\uparrow)$ and $(\downarrow\downarrow)$ Cooper pairs, respectively. The $z_\sigma$ term arises from Zeeman splitting in a magnetic field with $z \equiv z_\uparrow = -z_\downarrow$, the $c$ and $d$ terms  describe coupling between the two types of Cooper pairs, and the last term of Eq.~(\ref{eq:FreeEnergy}) accounts for the magnetic screening effect of the charge supercurrent. We point out that the standard form of the Ginzburg-Landau free energy~\cite{machida_incomplete_1978} only contains terms in which the total number of order parameters $\psi_\sigma$ and spatial derivatives $\grad$ does not exceed $4$. In this work, we include an additional symmetry-allowed term proportional to $d$ with $4$ order parameters and $2$ spatial derivatives that penalizes the charge supercurrent but not the spin supercurrent. As we will later find, this additional term is crucial for stabilizing HQVs. We also note that the dimensionless energy parameter $F_0$ is chosen such that $T_c F_0$ is the condensation energy of the entire superconductor with volume $V_0$ at zero temperature in the absence of a magnetic field ($\vec{a} = 0$) and any coupling ($c=d=0$).

\begin{figure}
\centering
\includegraphics[scale = 0.35]{Schematics.pdf}
\caption{Schematics of the Little-Parks experiment. The resistance of a thin superconducting ring with radius $R_0$, width $W \ll R_0$, and height $L$ is measured as a function of the applied magnetic field $\vec{H}$ near the superconducting critical temperature. The resistance itself is determined by applying a bias current $I_{\mathrm{bias}}$ and measuring the resulting voltage $\langle V \rangle$.}
\label{fig:Schematics}
\end{figure}

In the rest of this work, we focus on a thin superconducting ring of radius $R_0$ and width $W \ll R_0$ in a perpendicularly applied magnetic field $H$ (schematics shown in Fig.~\ref{fig:Schematics}). For $W \ll \xi$, where $\xi = \xi_0 (1-t)^{-1/2}$ is the finite-temperature coherence length, the order parameters $\psi_\sigma$ only depend on the polar angle $\theta$. Moreover, if $W$ is much smaller than the penetration depth $\lambda$, the screening effect of the charge supercurrent is negligible, and the vector potential in the symmetric gauge is simply given by $\vec{A} = \frac{1}{2} R_0 H \hat{\theta}$. Expressing the order parameters as $\psi_\sigma(\theta) = f_\sigma(\theta)e^{i\phi_\sigma(\theta)}$ in terms of the amplitudes $f_\sigma(\theta)$ and phases $\phi_\sigma(\theta)$, the Ginzburg-Landau equations derived from the free-energy functional of Eq.~(\ref{eq:FreeEnergy}) are then
\begin{widetext}
\begin{subequations}\label{eq:GL_LP}
\begin{eqnarray}
&& -\left[(1-t)+\frac{b_\sigma h}{r^2_0}\right]f_\sigma+f^3_\sigma+cf_\sigma f^2_{-\sigma}  = \frac{1}{r^2_0} \left\{ \frac{\partial^2f_\sigma}{\partial\theta^2}-f_\sigma\left[\frac{\partial\phi_\sigma}{\partial\theta}-h\right]^2 \nonumber \right\} \\
&& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \quad +\frac{d}{r^2_0} \left\{ f_{-\sigma}\frac{\partial^2(f_\sigma f_{-\sigma})}{\partial\theta^2}-f_\sigma f_{-\sigma}^2\left[\frac{\partial(\phi_\sigma+\phi_{-\sigma})}{\partial\theta}-2h\right]^2 \right\}, \label{eq:f_up}\\
&& \frac{\partial}{\partial\theta}\left\{f^2_\sigma\left[\frac{\partial\phi_\sigma}{\partial\theta}-h\right] + d(f_\sigma f_{-\sigma})^2\left[\frac{\partial(\phi_\sigma+\phi_{-\sigma})}{\partial\theta}-2h\right]\right\}= 0. \label{eq:phi_up}
\end{eqnarray}
\end{subequations}
\end{widetext}
Here $r_0 = R_0/\xi_0$ is a dimensionless ring radius, $b_\sigma = 2z_\sigma/\xi_0^2$ is a dimensionless Zeeman splitting, $h = H R_0^2\pi/\Phi_0$ is the number of flux quanta going through the ring, while $-\sigma$ indicates the opposite type of Cooper pair with respect to $\sigma$. We note that the conserved quantity within the curly brackets of Eq.~(\ref{eq:phi_up}) is a supercurrent that corresponds to the given type of Cooper pair [$(\uparrow\uparrow)$ or $(\downarrow\downarrow)$]; the charge and spin supercurrents are symmetric and antisymmetric combinations of these individual supercurrents, respectively. Since the order parameters $\psi_\sigma(\theta)$ must be single valued, the Ginzburg-Landau equations in Eq.~(\ref{eq:GL_LP}) are also supplemented with the boundary conditions $f_\sigma(\pi) = f_\sigma(-\pi)$ and $\phi_\sigma(\pi)-\phi_\sigma(-\pi) = 2\pi n_\sigma$, where $n_\sigma$ are arbitrary integers.

\section{Field-induced oscillations of the critical temperature}

In this section, we study the ``conventional'' Little-Parks oscillations~\cite{little_observation_1962} in the critical temperature of a spin-triplet superconducting ring as a function of the applied magnetic field. To describe these Little-Parks oscillations, we must establish the ground-state phase diagram of the system by enumerating and comparing ground-state candidates: stable solutions of Eq.~(\ref{eq:GL_LP}) that correspond to local minima of the free-energy functional in Eq.~(\ref{eq:FreeEnergy}). Such stable solutions take the general form of $f_\sigma (\theta) = \mathrm{constant}$ and $\phi_\sigma (\theta) = \phi^{(0)}_\sigma + n_\sigma \theta$, where $n_\sigma \in \mathbb{Z}$ are the fluxoid numbers for the two types of Cooper pairs and $\phi^{(0)}_\sigma$ are arbitrary reference phases~\cite{mccumber_time_1970}. Substituting this form into Eq.~(\ref{eq:f_up}), the constant values of the order parameters $f_\sigma$ are then solutions of the algebraic equations
\begin{eqnarray}\label{eq:solution_1stkind}
f_\uparrow \left[ f^2_\uparrow + \tilde{c} (n_\uparrow,n_\downarrow) f^2_\downarrow - \tilde{\alpha}_\uparrow (n_\uparrow) \right] &=& 0, \nonumber \\
\qquad f_\downarrow \left[ f^2_\downarrow + \tilde{c} (n_\uparrow,n_\downarrow) f^2_\uparrow - \tilde{\alpha}_\downarrow (n_\downarrow) \right] &=& 0,
\end{eqnarray}
where $\tilde{\alpha}_\sigma (n_\sigma) = (1-t)-r^{-2}_0[(n_\sigma-h)^2 - b_\sigma h]$ and $\tilde{c} (n_\uparrow,n_\downarrow) = c+r^{-2}_0 d (n_\uparrow+n_\downarrow -2h)^2$. These equations have four families of solutions that are listed in Table~\ref{tab:solutions} along with their free energies and physical interpretations. The trivial solution with $f_\uparrow = f_\downarrow = 0$ corresponds to a normal state, while the three nontrivial families describe superconducting states. For the $(n_\uparrow)_\uparrow$ and $(n_\downarrow)_\downarrow$ solutions, only one spin species [either $\uparrow$ or $\downarrow$] forms Cooper pairs with the other spin species remaining in a normal state. These solutions correspond to an effectively spinless superconductor that can only support full quantum vortices (FQVs) with a single integer fluxoid number ($n_\uparrow$ or $n_\downarrow$). In contrast, for the $(n_\uparrow,n_\downarrow)$ solutions, both spin species form Cooper pairs, and the result is a spin-triplet superconductor with both $(\uparrow\uparrow)$ and $(\downarrow\downarrow)$ pairing. Introducing the charge $n_c = (n_\uparrow + n_\downarrow)/2$ and spin $n_s = (n_\uparrow - n_\downarrow)/2$ fluxoid numbers, and recognizing that $n_c$ is the ``usual'' fluxoid number connected to the magnetic field, it is then clear that such a spin-triplet superconductor can support both FQVs (integer $n_{c,s}$) and HQVs (half-integer $n_{c,s}$).

\begin{table*}[t]
    \centering
    \renewcommand{\arraystretch}{2}
    \begin{tabular}{|c|c|c|c|c|}
        \hline Solution & $f_\uparrow^2$ & $f_\downarrow^2$ & Free energy ($F$) & Physical interpretation \\
        \hline Trivial & 0 & 0 & 0 & Normal (non-superconducting) state \\
        \hline $(n_\uparrow)_\uparrow$ & $\tilde{\alpha}_\uparrow$ & 0 & $-\frac{1}{2} F_0 \tilde{\alpha}_\uparrow^2$ & ``Spinless'' superconducting state with only ($\uparrow\uparrow$) Cooper pairs \\
        \hline $(n_\downarrow)_\downarrow$ & 0 & $\tilde{\alpha}_\downarrow$ & $-\frac{1}{2} F_0 \tilde{\alpha}_\downarrow^2$ & ``Spinless'' superconducting state with only ($\downarrow\downarrow$) Cooper pairs \\
        \hline $(n_\uparrow, n_\downarrow)$ & $\dfrac{\tilde{\alpha}_\uparrow - \tilde{c} \tilde{\alpha}_\downarrow} {1-\tilde{c}^2}$ & $\dfrac{\tilde{\alpha}_\downarrow - \tilde{c} \tilde{\alpha}_\uparrow} {1-\tilde{c}^2}$ & $-\frac{1}{2} F_0 \dfrac{\tilde{\alpha}_\uparrow^2 + \tilde{\alpha}_\downarrow^2 - 2\tilde{c} \tilde{\alpha}_\uparrow \tilde{\alpha}_\downarrow} {1-\tilde{c}^2}$ & Triplet superconducting state with ($\uparrow\uparrow$) and ($\downarrow\downarrow$) Cooper pairs \\ \hline
    \end{tabular}
    \caption{Four families of solutions for the order parameters $f_\sigma$ along with their respective free energies and physical interpretations. Each non-trivial family contains infinitely many distinct solutions labeled by the fluxoid numbers $n_\sigma$; the solutions for $f_\sigma$ and the corresponding free energies depend on $n_\sigma$ via $\tilde{\alpha}_\sigma (n_\sigma)$ and $\tilde{c} (n_\uparrow, n_\downarrow)$.}
    \label{tab:solutions}
\end{table*}

\begin{figure*}
\centering
\includegraphics[scale = 0.35]{Figure1.pdf}
\caption{Phase diagram of the ground-state fluxoid state (a) as a function of magnetic flux $h$ and coupling coefficient $d$ without Zeeman splitting ($b = 0$), (b) as a function of magnetic flux $h$ and temperature $t$ without Zeeman splitting, and (c) as a function of $h$ and $t$ with Zeeman splitting ($b \neq 0$). The white area marks the normal state, the yellow area marks a single-spin superconducting state with $f^2_\downarrow = 0$, and the red area marks a spin-triplet superconducting state with an integer fluxoid (FQV) around the ring ($n_{c,s} \in \mathbb{Z}$). In (a) \& (b), the light blue area indicates a degeneracy between two half-integer fluxoid (HQV) states with $n_s = \pm 1/2$, while in (c), the green (blue) color means that the HQV state with $n_s = +1/2$ ($n_s = -1/2$) is the ground state. In (c), the dashed and dash-dotted lines mark the upper and lower critical temperatures [Eqs.~(\ref{eq:tc_upper}) \& (\ref{eq:tc_lower})], respectively.}
\label{Figure1}
\end{figure*}

To determine the ground state of the ring, we need to compare the free energies of all solutions in Table~\ref{tab:solutions} while keeping in mind that each solution is only physical if $f_\uparrow^2 \geq 0$ and$f_\downarrow^2 \geq 0$. The resulting phase diagrams as a function of the temperature $t$, the magnetic field $h$, and the coupling constant $d$ are plotted in Fig.~\ref{Figure1}. We start understanding these results by comparing the various solutions within each family of Table~\ref{tab:solutions}. For the single-spin superconducting states $(n_\uparrow)_\uparrow$ and $(n_\downarrow)_\downarrow$ [marked by yellow color in Fig.~\ref{Figure1}], the free energy takes the exact form
\begin{equation}\label{eq:single-spin}
F_{(n_\sigma)_\sigma} = -\frac{F_0}{2} \left\{ (1-t) - \frac{1}{r^2_0} \left[ (n_\sigma - h)^2 - b_\sigma h \right] \right\}^2.
\end{equation}
We note that the expression inside the curly brackets must be positive for the given state to be valid. Hence, assuming $n \leq h \leq n+1$ (where $n$ is a non-negative integer) without loss of generality, the single-spin state with the lowest free energy at finite Zeeman splitting $b > 0$ (with $b \equiv b_\uparrow = -b_\downarrow$) is $(n)_\uparrow$ for $h < n+1/2$ and $(n+1)_\uparrow$ for $h > n+1/2$. For the spin-triplet states $(n_\uparrow,n_\downarrow)$, we expand the free energy up to $O(r_0^{-2})$ and express $n_\sigma$ in terms of $n_{c,s}$ to obtain
\begin{widetext}
\begin{equation}
F_{(n_\uparrow,n_\downarrow)} = F_0 \biggl\{-\frac{(1-t)^2}{1+c}+\frac{2(1-t)}{r_0^2}\frac{(1+c)n^2_s+[(1+c)+2(1-t)d](n_c-h)^2}{(1+c)^2}+O(r_0^{-4})\biggl\}. \label{eq:FreeEnergy_minima}
\end{equation}
\end{widetext}

Therefore, in the $r_0 \gg 1$ limit, we can restrict our attention to the FQV states with $n_s = 0$ and the HQV states with $n_s = \pm 1/2$ as all other states are penalized by the term $\propto n_s^2$. If we assume $n \leq h \leq n+1$ (with $n \in \mathbb{Z}$ and $n \geq 0$) again, the lowest-energy FQV state has $n_c = n$ for $h < n+1/2$ and $n_c = n+1$ for $h > n+1/2$. These two states correspond to $(n,n)$ and $(n+1,n+1)$ in the notation of Table~\ref{tab:solutions} and are both denoted by red color in Fig.~\ref{Figure1}. The two lowest-energy HQV states with $n_c = n+1/2$ and $n_s = \pm 1/2$, corresponding to $(n+1,n)$ and $(n,n+1)$, are degenerate up to $O(r_0^{-2})$. This degeneracy is split by a higher-order term in the free energy, $4 r_0^{-4} b h n_s (n_c - h) / (1 - c)$, such that the lowest-energy HQV state is $(n,n+1)$ for $h < n+1/2$ [blue color in Fig.~\ref{Figure1}] and $(n+1,n)$ for $h > n+1/2$ [green color in Fig.~\ref{Figure1}]. To compare the lowest-energy FQV and HQV states with one another, we finally recognize that the terms $\propto n_s^2$ and $\propto (n_c-h)^2$ have relative coefficients $(1+c)$ and  $(1+c)+2(1-t)d$ in Eq.~(\ref{eq:FreeEnergy_minima}). Given that $n_s^2 = 1/4$ for the HQV state and $0 \leq (n_c-h)^2 \leq 1/4$ for both states, the HQV state can only have lower energy than the FQV state if the coupling constant $d$ is positive. Thus, as previously stated, this coupling constant is crucial for stabilizing HQV states. Specifically, we find that the lowest-energy spin-triplet state is a HQV state within the field range $n+1/2-\Delta h/2 < h < n+1/2+\Delta h/2$ characterized by the width parameter
\begin{equation}
\Delta h = \frac{(1-t)d} {(1+c)+2(1-t)d} + O(r_0^{-2}). \label{eq:Delta-1}
\end{equation}
As shown by Eq.~(\ref{eq:Delta-1}) as well as Fig.~\ref{Figure1}, the field range in which a HQV state is energetically favorable increases as the temperature $t$ is lowered and as the coupling constant $d$ is increased. Conversely, this field range vanishes both at the critical temperature ($t \to 1$) and when the relevant coupling constant vanishes ($d \to 0$).

We are now ready to understand the field dependence of the critical temperature. We first recognize that, for $b h \neq 0$, there are in fact two critical temperatures; in addition to the upper critical temperature at which superconductivity first appears for one spin species, there is also a lower critical temperature marking a transition from a single-spin superconductor into a spin-triplet superconductor. These two critical temperatures can be determined by comparing the free energies of the lowest-energy single-spin and spin-triplet superconducting states [see Eqs.~(\ref{eq:single-spin}) and (\ref{eq:FreeEnergy_minima})] with each other and with the free energy of the normal state (which is zero). The upper critical temperature is found to be
\begin{equation}
t_{\mathrm{upper}} = 1 + \frac{1}{r^2_0} \left[ b|h| - (h - \lfloor h \rceil)^2 \right], \label{eq:tc_upper}
\end{equation}
where $\lfloor h \rceil$ is simply $h$ rounded to the nearest integer. This critical temperature has the same field-induced oscillation as the critical temperature of a spin-singlet superconductor~\cite{little_observation_1962} except for an additional linear increase with the field $|h|$ that shifts the maxima around integer $h$ [see Fig.~\ref{Figure1}(c)]. Up to $O(r_0^{-4})$, the lower critical temperature is given by
\begin{widetext}
\begin{equation}
t_{\mathrm{lower}} = 1 - \frac{[h - f(h)]^2 - c (h - \lfloor h \rceil)^2 + (1+c) b|h|} {r^2_0 (1-c)}- \frac{d [2h - \lfloor h \rceil - f(h)]^2 \{[h - f(h)]^2 - (h - \lfloor h \rceil)^2 + 2b|h|\}} {r^4_0 (1-c)^2} + O(r_0^{-6}) \label{eq:tc_lower}
\end{equation}
\end{widetext}
in terms of the piece-wise continuous function
\begin{equation}
f(h) = \bigg\{ \begin{matrix} \lfloor h \rceil \qquad \qquad \quad \mathrm{if} \,\, |h - \lfloor h \rceil| < (1 - \Delta h) / 2, \\ \lfloor 2h \rceil - \lfloor h \rceil \qquad \mathrm{if} \,\, |h - \lfloor h \rceil| > (1 - \Delta h) / 2, \end{matrix} \label{eq:fh}
\end{equation}
where the upper and the lower cases correspond to transitions into FQV and HQV states, respectively, and the field ranges around half-integer $h$ with transitions into HQV states are controlled by the width parameter
\begin{equation}
\Delta h = \frac{2db|h|} {r^2_0 (1-c)} + O(r_0^{-4}). \label{eq:Delta-2}
\end{equation}
In contrast to the upper critical temperature, the lower critical temperature shows an overall decrease with the field $|h|$ and contains additional maxima at all half-integer values of $h$ [see Fig.~\ref{Figure1}(c)] around which the transitions into the HQV states happen. We emphasize that, while Eqs.~(\ref{eq:Delta-1}) and (\ref{eq:Delta-2}) describe the same width parameter $\Delta h$, one cannot obtain Eq.~(\ref{eq:Delta-2}) by simply substituting Eq.~(\ref{eq:tc_lower}) into Eq.~(\ref{eq:Delta-1}) because the leading term of Eq.~(\ref{eq:Delta-2}) is $O(r_0^{-2})$ whereas Eq.~(\ref{eq:Delta-1}) is only accurate up to $O(1)$ terms. We further remark that a spin-singlet superconductor with multiple bands does not have two separate critical temperatures because its Ginzburg-Landau free energy includes bilinear coupling terms between different order parameters~\cite{erin_little-parks_2008}. Such coupling terms of the form $\propto \psi_{-\sigma}^* \psi_{\sigma}^{\phantom{*}}$ are forbidden in our case by spin-rotation symmetry around the field direction (i.e., the direction perpendicular to the ring).

\section{Magnetoresistance oscillations below the critical temperature}

In this section, we study the resistance of the superconducting ring in the spin-triplet state below the lower critical temperature in order to derive its oscillations as a function of the applied magnetic field. To do so, we will first consider thermal fluctuations in the order parameters $\psi_\sigma (\theta)$ that drive transitions between different free-energy minima. We will then use the computed transition rates to estimate the electrical resistance due to the resulting thermal decay of the charge supercurrent~\cite{mccumber_time_1970, cooper_resistance_2010}.

\begin{figure*}
\centering
\includegraphics[scale = 0.3]{Figure2.pdf}
\caption{(a,b) Saddle-point solutions for (a) the amplitude and (b) the phase of the superconducting order parameters $\psi_\uparrow$ and $\psi_\downarrow$ as the fluxoid numbers $(n_\uparrow,n_\downarrow)$ change from $(0,1)$  to $(1,1)$ at magnetic flux $h = 0.5$. (c) Free energies of the relevant stable states (solid lines) and the saddle points connecting them (dashed lines) as a function of the magnetic flux $h$.}
\label{Figure2}
\end{figure*}

During a thermal transition from one free-energy minimum to another one, the system goes through an appropriate free-energy saddle point, and the free-energy barrier controlling the transition rate is simply the free-energy difference between the saddle point and the original minimum. The two free-energy minima connected by the thermal transition correspond to two stable solutions of Eq.~(\ref{eq:GL_LP}) that differ in a fluxoid number $n_\sigma$. Therefore, the saddle-point solution must exhibit a phase slip in the corresponding order parameter $\psi_\sigma (\theta)$, i.e., a suppression of the amplitude $f_\sigma (\theta)$ around some angle $\theta$. For a spin-singlet superconductor, an analytical expression for this saddle-point solution was found in Ref.~\onlinecite{mccumber_time_1970}. Since the analytical solution only applies for a spin-triplet superconductor in the absence of coupling terms ($c=d=0$), we choose to solve Eq.~(\ref{eq:GL_LP}) numerically with a combination of a shooting method and an adaptive step size Runge-Kutta integration scheme.

An example of the numerically obtained saddle-point solution is shown in Figs.~\ref{Figure2}(a,b). This solution corresponds to a transition between two stable solutions with respective fluxoid numbers $(n_\uparrow,n_\downarrow) = (0,1)$ and $(1,1)$. Since $n_\uparrow$ changes from 0 to 1 while $n_\downarrow$ remains the same, $f_\uparrow$ goes through a strong suppression in a region of size $\xi = \xi_0 (1-t)^{-1/2}$ centered at $\theta=0$, whereas $f_\downarrow$ is only slightly perturbed around its stable constant solution.

Once the saddle-point configurations of $\psi_\uparrow$ and $\psi_\downarrow$ are found for a transition between $(n_\uparrow,n_\downarrow)$ and $(n'_\uparrow,n'_\downarrow)$, the saddle-point free energy associated with this transition, $F_{(n_\uparrow,n_\downarrow) \leftrightarrow (n'_\uparrow,n'_\downarrow)}$, can be calculated through Eq.~(\ref{eq:FreeEnergy}). For the stable solutions with fluxoid numbers $0 \leq n_{\sigma} \leq 1$ and the saddle-point solutions connecting them, the free energies are plotted in Fig.~\ref{Figure2}(c) as a function of the field $h$. We note that there is an exact degeneracy between $(n_\uparrow,n_\downarrow) = (0,1)$ and $(1,0)$ because we neglect the Zeeman splitting by setting $b=0$ in this section.

Since the free-energy barrier is simply given by $\Delta F_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)} = F_{(n_\uparrow,n_\downarrow)\leftrightarrow (n'_\uparrow,n'_\downarrow)}-F_{(n_\uparrow,n_\downarrow)}$ for the transition from $(n_\uparrow,n_\downarrow)$ to $(n'_\uparrow,n'_\downarrow)$, the thermally-activated rate of this transition can be estimated as~\cite{cooper_resistance_2010, halasz_fractional_2021}
\begin{equation}\label{eq:TransitionRate}
\Gamma_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)} \propto P_{(n_\uparrow,n_\downarrow)}\text{exp}\left[{-\beta\Delta F_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)}}\right],
\end{equation}
where $\beta = 1/T$ is the inverse temperature, and $P_{(n_\uparrow,n_\downarrow)}$ is the probability of the superconducting ring to be in the stable state $(n_\uparrow,n_\downarrow)$. At any finite temperature, this probability is given by
\begin{eqnarray}
P_{(n_\uparrow,n_\downarrow)} &=& \frac{1}{Z}\text{exp}\left[-\beta F_{(n_\uparrow,n_\downarrow)}\right],\nonumber \\
Z &=& \sum_{n_\uparrow,n_\downarrow}\text{exp}\left[-\beta F_{(n_\uparrow,n_\downarrow)}\right].
\end{eqnarray}
In the presence of a small bias current $I_{\mathrm{bias}}$ applied to a section of the ring, the free-energy barrier takes the modified form~\cite{cooper_resistance_2010} $\Delta \tilde{F}_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)} = \Delta F_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)} - (\delta n_\uparrow + \delta n_\downarrow) \Phi_0 I_{\mathrm{bias}}/4$ with $\delta n_\sigma \equiv n'_\sigma - n_\sigma$. The mean voltage $\langle V \rangle$ between the two end points of the section is then proportional to
\begin{widetext}
\begin{equation}
\langle V \rangle \propto \sum_{n_\uparrow,n_\downarrow} P_{(n_\uparrow,n_\downarrow)} \sum_{n'_\uparrow,n'_\downarrow} \left( \delta n_\uparrow + \delta n_\downarrow \right) \text{exp} \left[ -\beta \Delta F_{(n_\uparrow,n_\downarrow)\rightarrow (n'_\uparrow,n'_\downarrow)} + \beta \left( \delta n_\uparrow + \delta n_\downarrow \right) \frac{ \Phi_0 I_{\mathrm{bias}}} {4} \right].
\end{equation}
\end{widetext}
Finally, if we assume $I_{\mathrm{bias}} \ll T / \Phi_0$, the effective resistance takes the form
\begin{widetext}
\begin{equation}\label{eq:Resistance}
R_{\mathrm{eff}} = \frac{\langle V \rangle} {I_{\mathrm{bias}}} \propto \frac{\sum\limits_{n_\uparrow,n_\downarrow}\left\{\sum\limits_{\delta n_\uparrow = \pm 1}\delta n^2_\uparrow \text{exp} \left[ -\beta F_{(n_\uparrow,n_\downarrow) \leftrightarrow (n_\uparrow + \delta n_\uparrow, n_\downarrow)} \right] + \sum\limits_{\delta n_\downarrow = \pm 1}\delta n^2_\downarrow \text{exp} \left[ -\beta F_{(n_\uparrow,n_\downarrow) \leftrightarrow (n_\uparrow, n_\downarrow + \delta n_\downarrow)} \right] \right\}} {\sum\limits_{n_\uparrow,n_\downarrow} \text{exp} \left[ -\beta F_{(n_\uparrow,n_\downarrow)} \right]}.
\end{equation}
\end{widetext}
Here we ignore transitions that involve both fluxoid numbers $n_\sigma$ or change either fluxoid number by more than $1$ as they have larger free-energy barriers and their contributions to the resistance are negligible. We also disregard overall temperature-dependent prefactors from both the transition rate~\cite{mccumber_time_1970} and the expansion of $\langle V \rangle$ up to first order in $I_{\mathrm{bias}}$ as they only depend on the temperature via power laws, i.e., much weaker than the exponentials in Eq.~(\ref{eq:Resistance}).

\begin{figure*}
\centering
\includegraphics[scale = 0.3]{Figure3.pdf}
\caption{Magnetoresistance oscillations of a spin-triplet superconducting ring for different values of (a) the dimensionless energy parameter (i.e., condensation energy) $F_0$, (b) the coupling coefficient $d$, and (c) the dimensionless temperature $t$. The following parameter values apply to all three subfigures: $r_0 = 20$, $c = 0.1$, and $b = 0$. The resistance $R_{\mathrm{eff}}$ is calculated from Eq.~(\ref{eq:Resistance}) and normalized by its maximum value, $R_{\mathrm{eff}}^{\mathrm{max}}$, which is specified for each resistance curve.}
\label{Figure3}
\end{figure*}

Figure~\ref{Figure3} shows the magnetoresistance calculated from Eq.~(\ref{eq:Resistance}) at different values of the temperature $t$, the coupling constant $d$, and the energy parameter $F_0$. We readily observe the expected two-peak structure in the magnetoresistance oscillations~\cite{vakaryuk_effect_2011, yasui_little-parks_2017, cai_magnetoresistance_2022} and recognize clear trends in its behavior as a function of the three parameters. We first consider the separation between the two peaks that depends on $t$ and $d$ but not on $F_0$. As understood in previous works~\cite{vakaryuk_effect_2011, yasui_little-parks_2017, cai_magnetoresistance_2022}, each peak corresponds to a transition between a FQV state and a HQV state, with a HQV state being energetically favored in the narrow range between the two peaks. Hence, we immediately identify the separation between the two peaks as the width parameter $\Delta h$ in Eq.~(\ref{eq:Delta-1}) and understand that, for $d (1-t) \ll 1$, it is linearly proportional to both the coupling constant $d$ and the temperature difference $1-t$ with respect to the critical temperature.

We next focus on the width of each peak which, in conjunction with the peak separation, determines how well the two peaks are distinguishable from each other. In general, we observe that the peaks broaden if $F_0$ is decreased or $t$ is increased. To understand this result, we first notice from Eq.~(\ref{eq:Resistance}) that each peak extends over a field range in which the free-energy difference between the two lowest-energy fluxoid states $(n_\uparrow,n_\downarrow)$ is smaller than the temperature $T$; the resistance has a peak within such a field range because thermal phase slips between the two highly occupied lowest-energy states occur with an enhanced rate. Then, using Eq.~(\ref{eq:FreeEnergy_minima}) and assuming $r_0 \gg 1$ as well as $|1-t| \ll 1$, we find $|\partial F_{(n_\uparrow,n_\downarrow)} / \partial h| \sim F_0 (1-t) r_0^{-2}$ and estimate the peak width as $\delta h \sim T / |\partial F_{(n_\uparrow,n_\downarrow)} / \partial h| \sim r_0^2 T_c / [F_0 (1-t)]$. The conclusion is that the peak width is inversely proportional to both the energy parameter $F_0$ and the temperature difference $1-t$ with respect to the critical temperature. We note that the energy parameter $F_0$ is proportional to the total volume of the superconducting ring and can, in principle, be made arbitrarily large by increasing the out-of-plane ring height $L$; see the next section for quantitative estimates.

Even if the peaks are well separated and sufficiently narrow, the two-peak structure may still not be observable if the overall magnitude of the resistance is prohibitively small. According to Fig.~\ref{Figure3}, this problem occurs when $F_0$ is too large and/or $t$ is too small. From Eq.~(\ref{eq:Resistance}), the overall magnitude of the resistance is determined by the free-energy difference between the relevant stable and saddle-point solutions (i.e., the free-energy barrier). Since the saddle-point solutions have a localized suppression of superconductivity in a region of size $\xi = \xi_0 (1-t)^{-1/2}$ with respect to the stable solutions, this free-energy difference can be estimated as $\delta F \sim |F_{(n_\uparrow,n_\downarrow)}| \xi / R_0 \sim F_0 (1-t)^{3/2} r_0^{-1}$ if we assume $r_0 \gg 1$ and $|1-t| \ll 1$ again. The overall magnitude of the resistance is then found to be prohibitively small when $\beta \delta F \sim F_0 (1-t)^{3/2} / (r_0 T_c) \sim r_0 (1-t)^{1/2} / \delta h \gg 1$. Interestingly, this result reveals a fundamental trade-off between increasing the magnitude of the resistance and decreasing the peak width, which are both important for observing the two-peak structure. Still, Fig.~\ref{Figure3} demonstrates that a discernible two-peak structure with a reasonable magnitude of the resistance is possible with the right parameter values.

\section{Discussion and outlook}

In this work, we employed the Ginzburg-Landau approach to study the Little-Parks magnetoresistance oscillations of spin-triplet superconductors in the presence of HQVs harboring Majorana bound states. Focusing on a ring geometry with a sufficiently small width ($W \ll \xi, \lambda$), we first constructed the appropriate Ginzburg-Landau free energy and identified a specific higher-order term [the one proportional to $d$ in Eq.~(\ref{eq:FreeEnergy})] that is critical for stabilizing HQVs. Then, we used the associated Ginzburg-Landau equations to derive both the ``conventional'' Little-Parks oscillations of the superconducting critical temperature~\cite{little_observation_1962} and the closely related oscillations in the residual resistance due to thermal vortex tunneling below the critical temperature~\cite{cai_magnetoresistance_2022, sochnikov_large_2010, sochnikov_oscillatory_2010, cai_unconventional_2013, mills_vortex_2015}.

Our main result is a rigorous theoretical underpinning of the characteristic two-peak structure that is expected in the magnetoresistance oscillations~\cite{vakaryuk_effect_2011, yasui_little-parks_2017, cai_magnetoresistance_2022}. The two peaks demarcate a narrow field range of width $\Delta h$ [see Eq.~(\ref{eq:Delta-1})] in which it is energetically favorable to bind a HQV to the central hole of the superconducting ring. Since HQVs are stabilized by a higher-order term in the Ginzburg-Landau free energy, the peak separation $\Delta h$ is linearly proportional not only to the appropriate coupling constant $d$ but also to the relative temperature $1-t$ with respect to the critical temperature. In other words, the two-peak structure is observable in the magnetoresistance oscillations below the critical temperature but not in the ``conventional'' Little-Parks oscillations of the critical temperature itself.

We also established the general conditions under which the two-peak structure in the magnetoresistance oscillations is experimentally discernible. In particular, we identified a fundamental trade-off between making the peaks sufficiently narrow (hence, distinguishable from each other) and ensuring that the overall magnitude of the resistance is not overly suppressed. For the parameter values $r_0 = R_0 / \xi_0 \sim 10$ and $t = T / T_c \sim 0.9$ used in Fig.~\ref{Figure3}, we find that this trade-off corresponds to energy parameter $F_0 \sim 10^5$. Since $T_c F_0$ is defined as a total condensation energy at zero temperature, this quantity is on the order of $\nu \Delta^2 \sim \nu T_c^2 $, where $\nu$ is the electronic density of states at the Fermi level, and $\Delta$ is the superconducting pairing gap. Thus, we readily obtain $F_0 \sim \nu T_c \sim V_0 T_c / (a^3 E_F)$, where $V_0$ is the ring volume, $a$ is the lattice constant, and $E_F$ is the Fermi energy. For $T_c \sim 1$ K, $a \sim 1$ nm, and $E_F \sim 0.1$ eV, the ideal energy parameter $F_0 \sim 10^5$ then corresponds to ring volume $V_0 \sim 0.1$ $\mu$m$^3$, which is achieved, for example, by setting the ring radius, the ring width, and the ring height to $R_0 \sim 1$ $\mu$m, $W \sim 100$ nm, and $L \sim 100$ nm, respectively.

While the critical temperature marking the onset of superconductivity (i.e., the upper critical temperature) does not have a distinctive structure in its Little-Parks oscillations, we also found a lower critical temperature separating a single-spin superconducting state above and a spin-triplet superconducting state below. In turn, this lower critical temperature has two distinct minima around each half-integer value of $h$ [see Fig.~\ref{Figure1}(c)], which correspond to the two peaks in the related magnetoresistance oscillations~\footnote{Since the resistance decreases if the critical temperature is increased, it is not surprising that a two-peak structure in the resistance translates into a \emph{two-valley} structure in the critical temperature.}. To experimentally determine the lower critical temperature, one would need to selectively measure the electrical resistance of the spin species that is normal above and superconducting below the transition. For example, one could apply the bias current through half-metallic leads that can only emit or absorb one spin species but not the other one. We further note that this measurement would require a sharp superconducting transition with no residual resistance due to vortex tunneling below. As such, it would be in a different regime and need a much larger ring volume $V_0$ than discussed above. Since the radius and the width are restricted to $R_0 \lesssim 1$ $\mu$m and $W \lesssim 100$ nm by other considerations, this measurement would then correspond to a cylindrical geometry with a large height $L \gg 1$ $\mu$m as in the original Little-Parks experiment~\cite{little_observation_1962}.

Physically, the higher-order term proportional to $d$ in the Ginzburg-Landau free energy stabilizes HQVs by penalizing the charge supercurrent but not the spin supercurrent. In this sense, it plays an analogous role to the ratio of the spin and charge superfluid densities, $\gamma = \rho_{\mathrm{sp}} / \rho_{\mathrm{s}}$, in the London limit~\cite{chung_stability_2007, halasz_fractional_2021}. Indeed, from a comparison of Eq.~(\ref{eq:FreeEnergy_minima}) in this work and Eq.~(15) in Ref.~\onlinecite{halasz_fractional_2021}, we can identify the superfluid-density ratio as $\gamma = (1+c) / [(1+c) + 2(1-t)d]$, which is consistent with the general expectation that $\gamma \to 1$ at the critical temperature~\cite{vakaryuk_effect_2011}. Since it is also expected on general grounds that $\gamma < 1$ for interacting superconductors~\cite{chung_stability_2007, leggett_inequalities_1968, leggett_theoretical_1975}, we anticipate that, as a result of $d > 0$, our results apply to any spin-triplet superconductor in which $(\uparrow\uparrow)$ and $(\downarrow\downarrow)$ Cooper pairs are energetically favored over $(\uparrow\downarrow) + (\downarrow\uparrow)$ Cooper pairs. In the future, it would be interesting to understand how our results connect to the London limit and, in particular, how the fractional magnetoresistance oscillations in the presence of disorder predicted in Ref.~\onlinecite{halasz_fractional_2021} can be recovered from our general Ginzburg-Landau approach. This connection could then be used to establish the general conditions under which the fractional magnetoresistance oscillations are experimentally observable.

\section*{Acknowledgments}
This work was supported by the U. S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. This manuscript has been authored by employees of UT-Battelle, LLC under Contract No.~DE-AC05-00OR22725 with the U.S. Department of Energy. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

\end{document}


\section{Introduction}
\label{sec1}

Many geophysical and astrophysical flows are confined in thin
layers of small aspect ratio either by material boundaries or by
physical mechanisms which constrain the motion, such as rotation or
stratification.
The vertical extension (thickness) of such layers can be much
smaller than the typical horizontal scales, while being at the same
time much larger than the dissipative viscous scales. As a consequence,
turbulent flows in those quasi-two-dimensional geometries display a
rich phenomenology with both three-dimensional (3D) features at small scales
and two-dimensional (2D) properties at large scales.

Remarkably, numerical simulations have shown that a physical
confinement is not necessary to observe a two-dimensional phenomenology.
Even fully periodic simulations in a box with large aspect ratio
$L_z \ll L_x=L_y$, forced at intermediate scales,
produce a split energy cascade in which a
fraction of the energy flow to large scales (as in a pure 2D flow) and the
remaining part goes to small scales producing the 3D direct cascade
\cite{smith1996crossover,celani2010turbulence,alexakis2018cascades}.
In this case, the key parameter which controls the relative flux of
energy in the two cascades is the ratio $S=L_z/L_f$ between the confining
scale and the characteristic scale of the forcing $L_f$.
In the limit $S \to 0$ of very small thickness, vertical motion is
suppressed by viscosity and the flow fully recovers the 2D
phenomenology
\cite{musacchio2017split}.
The phenomenology of this energy split can be changed by other physical
ingredients, including solid body rotation \cite{deusebio2014dimensional}
and stratification \cite{sozza2015dimensional}. Moreover, the
effects of the transition are important for several statistics, including
the chaoticity of the flow \cite{clark2021chaotic}.
Since at larger scale the flow becomes more two-dimensional,
the inverse cascade of energy in thin layer proceeds until it reaches
the largest scale and produces a large-scale structure, called the
condensate \cite{xia2009spectrally,laurie2014universal,musacchio2019condensate,van2019condensates,fang2021spectral}.

The phenomenology changes in the presence of confining physical
boundaries, as in the case of laboratory experiments performed in a
thin layer of fluid confined by gravity
\cite{shats2010turbulence,byrne2011robust,xia2011upscale}.
The bottom (and lateral) boundary of the tank produces a boundary
layer which dissipates a relevant fraction of the energy injected in the system
\cite{van2009studies,shats2010turbulence} and thus reduces the turbulent
flux, in particular in the case of a single layer of
fluid \cite{byrne2011robust}.
Experiments with a double layer, in particular of immiscible fluids,
reduce the damping rate induced by the bottom wall, and produce an inverse
cascade of energy \cite{chen2006physical,byrne2011robust}.

In this work we systematically study the turbulent flow in a
thin layer by extensive direct numerical simulations
of the Navier-Stokes equations in a box with no-slip
boundary conditions (BC) at the bottom, and a free-slip
BC at the top, similarly to laboratory experiments of
free-surface turbulence.
Simulations are done at high resolution
and large Reynolds numbers which allows the development of a
fully 3D turbulent motion. The forcing scale $L_f$ is fixed, and
we vary the ratio $S$ by changing the thickness $L_z$ of the box.
We find that, in the range of parameters explored here, the thin
layer is unable to sustain an inverse cascade of energy.
We observe a transient inverse cascade, of duration which depends
on the thickness $S$, but the 3D motion in the bottom boundary layer
eventually propagates to the full layer and the flow becomes fully
three-dimensional with the inverse cascade being suppressed.
We also study the decaying regime of our systems, and we
find a complex behavior which cannot be described by a simple
exponential law.

The remaining of this paper is organized as follows. In Sec.~\ref{sec2},
we present the details of the numerical simulations. Section~\ref{sec3}
is devoted to the evolution of the global quantities of the flow,
while Sec.~\ref{sec4} discusses small scale statistics and in particular
the presence of a direct or inverse cascade. In Sec.~\ref{sec5}, we report
the results of the decaying simulations and finally Sec.~\ref{sec6} is
devoted to the conclusions.

\section{Models and direct numerical simulations}
\label{sec2}

We consider the 3D Navier-Stokes equations for an incompressible
velocity field ${\bm u}({\bm x},t)=(u,v,w)$
(with ${\bm \nabla} \cdot {\bm u}=0$) in a domain of dimension
$L_x \times L_y \times L_z$
\begin{equation}
{\partial {\bm u} \over \partial t} + {\bm u} \cdot {\bm \nabla} {\bm u} =
- {\bm \nabla} p + \nu \nabla^2 {\bm u} + {\bm f},
\label{eq1}
\end{equation}
where the constant density has been adsorbed into the pressure $p$ and
$\nu$ is the kinematic viscosity.
The two-dimensional forcing ${\bm f}$ is restricted to the two
horizontal ($x,y$) components (2D2C) ${\bm f}({\bm x})=(f_x(x,y),f_y(x,y),0)$.
It is Gaussian, white in time, and in Fourier space is confined
in a narrow cylindrical shell of wavenumbers centered around
$k_f=2 \pi/L_f=8$. One reason to have a $2D$ forcing is that it is
independent on the thickness of the flow and therefore we use
the same forcing for all the simulations at different $L_z$.
Moreover, thanks to the delta-correlation in time, the rates of injection of
energy $\varepsilon$ is fixed and does not depend on $L_z$ or
on the properties of the flow.

Boundary conditions are periodic in the horizontal direction $(x,y)$
while, to simulate free-surface turbulence, we impose a no-slip BC
at the bottom $z=0$ and a free-slip BC at the top
$z=L_z$. We have therefore $u=v=w=0$
at $z=0$
while $\partial_z u=\partial_z v=0$ and $w=0$ at $z=L_z$.
We remark that these BC have been previously used for
numerical studies of free-surface turbulence \cite{nagaosa1999direct}.

The equations of motion are solved numerically with a uniform spacing in all
directions. To solve the problem, we use the flow solver \textit{Fujin}, an
in-house code, extensively validated and used in a variety of
problems~\cite{olivieri2020dispersed, rosti2021shear, mazzino2021unraveling,
brizzolara2021fiber, mazzino2022puff, hori2022eulerian}, based on the
(second-order) finite-difference method for the spatial discretization and the
(second-order) Adams-Bashforth scheme for time marching. See also
\texttt{https://groups.oist.jp/cffu/code} for a list of validations.
Simulations are performed at a fixed horizontal resolution and varying
the vertical resolution depending on $S$ with a constant viscosity and
energy input (see Table~\ref{table1}).
Additional simulations at different viscosities (not
discussed here) produced similar results.
All the simulations start from an initial zero-velocity field
and reach a statistically stationary states characterized
by a constant energy. Another set of simulations starts from
these asymptotic states and study the decaying regime by switching
off the forcing.
We remark that since $\langle f_i \rangle=0$
(the average is defined over the $(x,y)$ planes) we have no mean flow
$\langle u_i \rangle=0$ and we will consider the statistics of
the fluctuating velocity field only.

\begin{table}[h!]
\begin{tabular}{cccccc}
$Run \#$ & $N_z$ & $S$ & $E_{2D}$ & $E_z$ & $Re$ \\ \hline
$1$ & $16$ & $0.03125$ & $3.9$ & $0.001$ & $1580$ \\
$2$ & $32$ & $0.0625$ & $7.8$ & $0.07$ & $2250$ \\
$3$ & $48$ & $0.09375$ & $9.9$ & $0.22$ & $2550$ \\
$4$ & $64$ & $0.125$ & $11.6$ & $0.41$ & $2780$ \\
$5$ & $96$ & $0.1875$ & $12.6$ & $0.68$ & $2920$ \\
$6$ & $128$ & $0.25$ & $13.9$ & $0.96$ & $3100$ \\
$7$ & $160$ & $0.3125$ & $12.4$ & $1.1$ & $2940$ \\
$8$ & $512$ & $1.0$ & $12.5$ & $2.1$ & $3060$  \\
$9$ & $1024$ & $2.0$ & $12.5$ & $2.7$ & $3100$
\end{tabular}
\caption{Parameters of the simulations.
$N_z$ is the resolution in the $z$ direction, $S=L_z/L_f$ the thickness in the $z$
direction (in units of the forcing scale $L_f$), $E_{2D}$ and $E_z$ are the
energies of the horizontal and vertical component of the velocity in the
stationary state (the total energy is $E=E_{2D}+E_z$). The Reynolds
number is defined as $Re=\sqrt{E} L_f/\nu$.
For all the runs $L_x=L_y=2 \pi$ with
a resolution $N_x=N_y=4096$. The viscosity is $\nu=9.8 \times 10^{-4}$
and the random forcing is active on the scale $L_f=L_x/8$ with a fixed
energy input $\varepsilon_{I}=50$.
}
\label{table1}
\end{table}

\section{Large-scale properties of the flow}
\label{sec3}

The $2D2C$ random forcing initially produces a two-dimensional flow.
Since it cannot transfer energy to small scale,
energy dissipation during this first phase is negligible and the
kinetic energy grows approximatively as $E(t) \simeq \varepsilon_I t$.
After this phase, vertical motions start to develop and eventually produce a
three-dimensional turbulence with a transfer of energy to small scales where
it is dissipated. As a consequence, the kinetic energy of the flow
reaches a (statistically) stationary state.
This is clearly shown in Fig.~\ref{fig1} where we plot the time
evolution of the total kinetic energy $E$ and of the vertical component
$E_z$ for the Run $6$.
For $t \le 0.3$, the vertical kinetic energy is negligible and the
total energy grows at the input rate. After $t \simeq 0.4$, vertical
motion sets in and the turbulent transfer to small scales produces a
viscous dissipation, which reduces the energy which eventually reaches
a stationary state. We remark that in this stationary state the vertical
kinetic energy is still much smaller than the horizontal one
(see Fig.~\ref{fig1} and Table~\ref{table1}).

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig1.eps}}
\caption{Time evolution of the total kinetic energy $E$ (red line) and
of the vertical kinetic energy $E_z$ (blue line) for the forced run $6$
at $S=0.25$. The black dashed line represents the growth rate due to
the energy input, $\varepsilon_I t$.
Inset: Time evolution of the total kinetic energy $E$ for all the runs.
Thickness $S$ (and Run number) increases from bottom (red) to top (black)
lines. The black dashed line represents $\varepsilon_I t$.}
\label{fig1}
\end{figure}

From Fig.~\ref{fig1} we can clearly distinguish two phases in the
turbulent flow: a two-dimensional regime at initial times and a
three-dimensional regime at late times. This picture is observed
(with quantitative differences) for all the simulations at different
thickness as shown in the inset of Fig.~\ref{fig1}.
The fact that kinetic energy reaches a constant value indicates that,
for any thickness, the flow is unable to sustain an inverse cascade (which
would keep the energy increasing). The reason why the {\it energy split}
scenario (in which both a direct and an inverse energy cascades are present)
is not observed in our simulations is the main result of our work and
will be discussed in details below.

From Fig.~\ref{fig1} we observe that the asymptotic value of the energy
grows with the thickness, more rapidly for lower values of $S$.
This is shown more clearly in Fig.~\ref{fig2}
together with the dependence of $E_z/E$ on $S$ in stationary
conditions. We see that for
$S \lesssim 0.3$ the asymptotic energy has a strong dependence on $S$,
while for larger values of the thickness, $E$ reaches an almost constant
plateau.
The inset of Fig.~\ref{fig2} displays the dependence of the ratio
$E_z/E$ with the thickness.
In this case we observe a growth for all the values of
$S$, indicating that the presence of the bottom
boundary affects the vertical motion even at $S \approx 1$.
We remark that the flow remains
anisotropic also at large $S$ since $E_z/E<1/3$.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig2.eps}}
\caption{Asymptotic kinetic energy $E$ as a function of the thickness $S$.
Inset: ratio $E_z/E$ as a function of $S$.
}
\label{fig2}
\end{figure}

Since the flow has no mean velocity,
$\langle {\bf u} \rangle = 0$, it is natural to consider the vertical
profiles of the rms velocities, i.e.
${\bm u}_{rms}(z)=\langle {\bm u}({\bf x})^2 \rangle^{1/2}$. The different
boundary conditions on the horizontal $(u,v)$ and vertical $w$ components,
together with the $2D$ forcing, produce different profiles for the
horizontal and vertical components.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig3.eps}}
\caption{Vertical profile of the horizontal velocity
fluctuations $u_{rms}$. The inset shows the behavior close to the
boundary at $z=0$. Colors as in the inset of Fig.~\ref{fig1}.}
\label{fig3}
\end{figure}

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig4.eps}}
\caption{Vertical profile of the vertical velocity
fluctuations $w_{rms}$. Colors as in the inset of Fig.~\ref{fig1}.}
\label{fig4}
\end{figure}

Fig.~\ref{fig3} shows the vertical profiles of one component of
the horizontal velocity $u_{rms}$ in stationary conditions.
We observe, close to $z=0$, a boundary layer region, strongly affected
by the presence of the bottom wall, where velocity fluctuations increase
rapidly with $z$. For larger values of $z$ (and not too small
$S$) horizontal velocity fluctuations saturates to an approximately
constant value which increases with $S$.
In these cases, we observe that close to the upper boundary $z=L_z$,
the horizontal velocity fluctuations increase. This is a consequence
of the free-slip boundary conditions and has been already observed in
free-surface channel flow \cite{nagaosa1999direct,lovecchio2013time}.
A possible explanation of this effect is obtained from a
Taylor expansion of the velocity field close to the surface
$z=L_z$. Introducing, for simplicity of notation,
 the shifted variable $Z=L_z-z$ (such that
the free surface is at $Z=0$), boundary conditions imply for
one horizontal component of the velocity
$u({\bf x})=u_0(x,y)+u_2(x,y)Z^2 + O(Z^4)$.
Therefore, the variance of the velocity close to the surface has the
expression
\begin{equation}
\langle u^2 \rangle=\langle u_0^2 \rangle + 2 \langle u_0 u_2 \rangle Z^2
+ O(Z^4).
\label{eq3.1}
\end{equation}
Assuming that the energy dissipation rate averaged over
horizontal planes is independent on $z$ (which is verified in our simulation
for values of $z$ not too close to the bottom plane), we expect that
$\langle \varepsilon \rangle = -a \nu \langle u \nabla^2 u \rangle > 0$,
where the positive constant $a$ depends on the details of the flow
($a=3$ for isotropic turbulence).
By Taylor expansion we can write
\begin{equation}
u \nabla^2 u \simeq (u_0 + u_2 Z^2)(\partial_x^2+\partial_y^2+\partial_Z^2)
(u_0+u_2 Z^2) \, .
\label{eq3.2}
\end{equation}
Using integration by part in the periodic directions $x$ and $y$ we obtain
\begin{equation}
\langle u \nabla^2 u \rangle = - \langle (\partial_x u)^2 \rangle -
\langle (\partial_y u)^2 \rangle + 2 \langle u_0 u_2 \rangle + O(Z^2)
\simeq - \varepsilon/\nu < 0.
\label{eq3.3}
\end{equation}
This expression suggests that, in the absence of cancellations of leading
terms, $\langle u_0 u_2 \rangle \propto - \varepsilon/\nu < 0$
and therefore, from (\ref{eq3.1}), that $\langle u^2 \rangle$
deceases moving away from the free surface as observed in
Fig.~\ref{fig3}.
A parabolic fit of the velocity variance close to $z=L_z$
is quantitatively consistent with the above predictions.
Close to the bottom boundary with no-slip BC we observe that the
extension of the boundary layer region is weakly dependent on $S$, as it is
shown in the inset of Fig.~\ref{fig3}.

Vertical profiles of the vertical velocity fluctuations are shown
in Fig.~\ref{fig4}. At variance with the horizontal components,
here the velocity (and therefore its fluctuations) vanishes also at
the upper free surface. The maximum value of fluctuations is observed
approximately in the middle of the domain, even if, due to the
different boundary conditions, the profiles are not symmetric with
respect the central plane $z=L_z/2$.
The ratio of the lines plotted in Fig.~\ref{fig4} and in Fig.~\ref{fig3}
gives the profile of the anisotropy of the velocity field. Even for the
largest values of $S$, we have that in the central part of the domain
$w_{rms}/u_{rms} \simeq 0.5$.

\section{Small-scale statistics and transient inverse cascade}
\label{sec4}

As discussed in the previous Section, the 2D2C forcing produces
initially a two-dimensional flow. This flow is non-stationary,
as shown in Figs.~\ref{fig1}, and eventually develops
instabilities which result in a three-dimensional motion.
In order to better understand and characterize this transition,
we study the small scale statistics of the turbulent flow at
different times and horizontal planes.
We consider here the run at $S=2$ for which the results are more clear,
but a similar scenario is observed for the other thickness also.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig5a.eps}
\includegraphics[scale=0.75]{fig5b.eps}}
\caption{2D energy spectra for the simulation at $S=2$ computed at
$z=L_z/4$ (left) and $z \simeq L_z$ (right) at different times,
from $t=0.2$ (light red)
to time $t=1.2$ (black) with steps $\Delta t=0.2$.
The dashed line represents the
Kolmogorov spectrum $k^{-5/3}$.}
\label{fig5}
\end{figure}

Figure~\ref{fig5} shows the horizontal spectra (computed on the
($x,y$) plane) of the
full velocity field in two different depth $z$ and at different times.
Both at $z=L_z/4$ and at $z \simeq L_z$ (i.e. close to the lower and
upper boundary respectively) the flow initially (at $t=0.2$)
develops fluctuations at small scales (i.e. at $k>k_f=8$) with a
steep spectrum compatible with a 2D direct cascade of enstrophy.
During this first
stage, in which the energy increases approximately with the
energy injection (see Fig.~\ref{fig1}), the flow transfers some
energy at large scale as prescribed by a two-dimensional phenomenology.
At $t=0.4$ the flow close to the bottom has already developed
a Kolmogorov scaling $k^{-5/3}$, compatible with a 3D direct cascade,
while the spectrum of the flow close to the surface is still steep
(second line from bottom in both plots).
At later times $t \ge 0.6$, spectra display a Kolmogorov scaling
in both planes.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig6.eps}}
\caption{2D energy spectra from the simulation at $S=2$ at time $t=0.4$
computed at different planes from $z=L_z/8$ (clear red) to
$L_z \simeq L_z$ (black).
The dashed line represents the Kolmogorov spectrum $k^{-5/3}$.}
\label{fig6}
\end{figure}

The interpretation of these results is that two-dimensional turbulence,
which is initially produced by the 2D2C forcing, is transient and the
flow develops a three-dimensional direct cascade starting from the layers
closer to the bottom boundary. In this sense, the
transition from 2D to 3D flow is not uniform in space, but
proceeds from the bottom layer towards the top layer. This is
clearly shown in Fig.~\ref{fig6}, where we plot the energy spectra
computed on different horizontal planes in the domain at the
intermediate time $t=0.4$. It is clear that, while the upper layers
(at $z \simeq L_z$) have a steep spectra compatible with a direct
2D cascade, the lower layers close to the bottom have already
developed a 3D cascade with a Kolmogorov spectrum.

The transition from 2D to 3D dynamics at different layers is confirmed
by the analysis of structure functions (SF) in physical space,
in particular by the third-order SF which contains
information about the flux of energy \cite{boffetta2012two}.

From the longitudinal velocity increments
$\delta u_{L}(\ell,{\bm x},t) \equiv
({\bm u}({\bm x}+{\bm \ell},t)-{\bm u}({\bm x},t)) \cdot
{\bm \ell}/\ell$
where ${\bm \ell}$ is a vector on the $(x,y)$ plane,
we define the SF of order $p$ as
\begin{equation}
S_p(\ell;z) = \langle \delta u_{L}^p(\ell,{\bm x},t) \rangle
\label{eq4.1}
\end{equation}
where, as in Section~\ref{sec1}, the average is over the plane
$(x,y)$ and time.
We remark that three-dimensional turbulence is characterized by a
negative third-order SF corresponding to a direct cascade of
turbulent fluctuations to small scales. In particular, in
3D homogeneous-isotropic turbulence one has
$S_3(\ell)=-(4/5) \varepsilon \ell$
\cite{frisch1995turbulence}.
In two dimensions one has, on the contrary, an inverse cascade of
turbulent fluctuations at scales larger than $L_f$
with a positive third-order SF given,
in homogeneous-isotropic conditions, by
$S_3(\ell)=(3/2) \varepsilon \ell$
\cite{boffetta2012two}.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig7.eps}}
\caption{Horizontal longitudinal third-order velocity structure functions
$S_3({\ell;z})$
computed at three times of the simulation, $t=0.2$ (red lines),
$t=0.4$ (blues lines) and $t=0.6$ (green lines).
Continuous lines correspond to $z=3/4 L_z$, dashed line to
$z=1/2 L_z$ and dashed-dotted line to $z=1/4 L_z$.
Simulation at $S=2$.}
\label{fig7}
\end{figure}

In Fig.~\ref{fig7} we plot the horizontal third-order longitudinal velocity
SF at three different times and at three different depth corresponding
to $z/L_z=1/4,1/2,3/4$. At short time, $t=0.2$ the SF is positive,
corresponding to an inverse energy cascade, at all the three depths considered,
consistent with the spectrum shown in Fig.~\ref{fig4}.
At the intermediate time $t=0.4$ (which still correspond to the
growing phase of the total energy, see inset of Fig.~\ref{fig1}), the SF in the
upper layer (continuous line) is still positive, while it becomes
negative in the lower layer close to the bottom boundary (dotted-dashed
line). At the intemediate layer the SF change sign with the scale and
at large scales is still positive.
This confirms the picture in Fourier space observed in Fig.~\ref{fig6}.
At time $t=0.6$, which corresponds to the peak of the kinetic
energy in Fig.~\ref{fig1}, the SFs become negative for all the
depth considered, indicating a complete transition to 3D turbulence
in the whole domain.
This is at variance to what observed in simulations in fully
periodic domains where, for intermediate values of $S$, it is observed an
inverse cascade at large scales together with a direct cascade at
small scales, both in stationary conditions \cite{celani2010turbulence}.
Remarkably, laboratory experiments in a conducting fluid
show a similar phenomenology with the transion from positive to
negative third-order SF \cite{gledzer2011structure}.

The physical interpretation of these results, which are qualitatively
confirmed also for the other simulations at smaller $S$, is that
the 2D2C forcing initially produces a quasi two-dimensional flow,
which is turbulent (positive $S_3(\ell)$) in the whole domain.
As the energy increases, the friction on the bottom boundary
induces vertical motions which make the flow three-dimensional
(negative $S_3(\ell)$) starting from the lower layers.
Eventually, at longer times, the whole flow becomes three dimensional,
turbulent fluctuations are dissipated by viscosity and the kinetic
energy decreases to reach a stationary state.

\section{Decaying turbulence}
\label{sec5}

In order to better understand the effects of the bottom layer on the
turbulent flow, we performed additional simulations in decaying conditions,
i.e. by integrating (\ref{eq1}) without the forcing term ${\bm f}$.
The initial conditions for the decaying simulations are taken from the
forced runs in stationary conditions, i.e. from the last time of
the inset of Fig.~\ref{fig1}.

\begin{figure}[h!]
\centerline{\includegraphics[scale=0.75]{fig8.eps}}
\caption{Decaying of total kinetic energy in unforced simulations at different
thickness $S=0.03125$ (red lower line), $S=0.0625$ (blue intermediate line),
$S=0.125$ and $S=0.3125$ (green and orange upper lines).
Inset: local exponential slope of $\log (E(t)/E(0))$ for simulations
at $S=0.03125$ (red upper line), $S=0.0625$ (blue intermediate line),
$S=0.125$ (green lower line). Dashed lines represent the viscous decay
rate $\lambda=\pi^2 \nu/(4 L_z^2)$.}
\label{fig8}
\end{figure}

The evolution of the energy (normalized with the initial energy $E(0)$) is
shown in Fig.~\ref{fig8} (where the initial time $t=0$ is now the time at
which the forcing is switched off). From a qualitative point of view, it is
evident that the decay is faster for thinner layers, indicating the
importance of the bottom boundary for the dissipation of energy.
Nonetheless, we see that energy dissipation rate becomes almost independent in
the case of thick layers (the lines for $S=0.125$ and
$S=0.3125$ are practically identical).

The lin-log plot of Fig.~\ref{fig8} suggests that, while in the thinnest case
$S=0.03125$ the long-time decay is with good approximation exponential,
this is not the case for the
simulations at larger values of $S$. This is quantitatively confirmed by the
inset of Fig.~\ref{fig8} where we plot the local rate of exponential decay for
three cases. The rate $2 \lambda$ is obtained by plotting the local slope of
$\log (E(t)/E(0))$ vs time, while the dashed lines represent the
theoretical viscous decaying rate of a 2D flow with friction, given by
$\lambda=\nu \pi^2/(4 L_z^2)$
\cite{shats2010turbulence}.
It is clear that, while the flow at $S=0.03125$, after a short transient,
reaches the exponential decay with the predicted viscous rate, the other cases
with larger $S$ display a more complex decay law which cannot be
simply described by an exponential law and a corresponding  friction
coefficient $\lambda$. A complex decaying law in quasi 2D experiments
has been recently reported and interpreted as different stages of
exponential decay with different decaying constant \cite{fang2017multiple}.
The results shown in the inset of Fig.~\ref{fig8} indicate that in our case
it is difficult to recognize a, even transient, exponential regime and
that the interplay of 2D and 3D motion produces a complex decaying
phenomenology.

\section{Conclusions}
\label{sec6}

We studied the dynamics and statistics of turbulence in a thin layer
with no-slip BC on the bottom surface, forced by a 2D2C forcing.
For the range of thickness explored, we find that the flow is
unable to sustain an inverse energy cascade: after a initial transient,
the flow develops a three-dimensional direct cascade which starts from the
bottom and eventually propagates in the whole domain.
The analysis of the energy spectrum and of the third-order structure
function shows that at intermediate times, a 2D like and a 3D like
phenomenologies can coexist at different depth in the flow.
This is in contrast to what observed in homogeneous simulations in
the absence of boundaries where the
the flow can sustain simultaneously an inverse cascade of energy to
large scales and a direct cascade to small scales \cite{celani2010turbulence}.
Moreover, our results are also at odds with several laboratory
experiments where an inverse cascade is observed
\cite{paret1997experimental,chen2006physical,von2011double,byrne2011robust,xia2011upscale},
but mostly in the presence of two layers of fluids (miscible
or not). The differences in these cases are probably due to the fact
that the dynamics in the upper layer (where the flow is studied), is
partially decoupled from the lower layer which
is affected by the no-slip boundary conditions.

Finally, we studied the decaying behavior of the thin turbulent layer.
Also in this case we find differences with respect to laboratory
experiments where the decay of kinetic energy is exponential and
therefore it is parameterized by a single friction coefficient
\cite{shats2010turbulence}.
We observe a clear exponential decay only
for the simulation with the thinnest layer, while in the other cases
a more complex decay law is observed.
The different behavior in this case can be not only ascribed to the
presence of stratification, but also to the
different initial flow conditions between experiments and simulations.

In this work we have changed only one parameter of the flow (the
thickness $S$), while there are other variables which can produce
a different phenomenology, mainly the viscosity (i.e. the Reynolds
number) and the forcing statistics \cite{poujol2020role}.
Further work is therefore needed to
uncover the reach phenomenology of turbulent thin layers.

\section*{Acknowledgments}
We acknowledge HPC CINECA for computing resources
(INFN-CINECA grant no. INFN22-FieldTurb).
M.E.R. is supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan. M.E.R. also acknowledges the computational time provided by the Scientific Computing section of Research Support Division at OIST.