\section{Fitting parameters of empirical function} \label{sec_app1}
As briefly explained in Sect.~\ref{sec_ah}, the behavior of the absolute humidity dependence of the calibration factor is well described by an empirical function, $C_{F}=A-B\sqrt{A_{H}}$. Table~\ref{table_fit} summarizes the fitting parameters obtained in this study as well as the past studies.

\begin{table*}
\centering
{\scriptsize
\caption{Summary of the fitting parameters of the empirical function. Those parameters are obtained by filling purified air.
\label{table_fit}
\begin{tabular}{cccccc
\hline
Detector & Size of  & High voltage & $A$ & $B$ & Comment \\
Volume & PIN-photodiode & [kV] & & & \\
\hline
$1$~L & $10\times10~\mathrm{mm}$ & $-0.12$ & $(12.86\pm0.40)\times10^{-4}$ & $(1.66\pm0.19)\times10^{-4}$ & Ref.~\cite{Pronost:2018ghn} \\
$70$~L & $16\times16~\mathrm{mm}$ & $-1.5$ & $2.03$ & $1.08$ & Ref.~\cite{TAKEUCHI1999334} \\
 & & & & & Constant~(0.73)  \\
  &  & & & & above $1.6~\mathrm{g/m^{3}}$ \\
$80$~L & $18\times18~\mathrm{mm}$ & $-2.0$ & $2.25\pm0.05$ & $0.29\pm0.05$ & Ref.~\cite{Nakano:2017rsy} \\
$80$~L & $18\times18~\mathrm{mm}$ & $-2.0$ & $2.39\pm0.04$ & $0.42\pm0.05$ & This study \\
$80$~L & $28\times28~\mathrm{mm}$ & $-2.0$ & $2.48\pm0.04$ & $0.37\pm0.05$ & This study \\
\hline
\end{tabular}
}
\end{table*

\section{Evaluation of the background rate including the calibration setup} \label{sec_bg_setup}

We measured the background rate by circulating the purified air in the calibration setup in order to evaluate some contributions of Rn emanation from the other components of the calibration setup.  Figure~\ref{fig_bg_spect2} shows the measured background spectrum when the filled purified air is circulated. The air flow rate was controlled at $0.6$~L/min by the mass flow controller.

\begin{figure*}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{hist_bg18_app.eps}
\end{minipage}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{hist_bg28_app.eps}
\end{minipage}
\caption{The daily background spectrum measured by two Rn detectors by circulating the purified air. The left~(right) panel shows the spectrum measured by the Rn detector equipped with $18\times18~\mathrm{mm}$ ~($28\times28~\mathrm{mm}$) PIN-photodiode.}
\label{fig_bg_spect2}
\end{figure*}

Figure~\ref{fig_bg_rate2} shows the daily count rate in this measurement. The measured background rate of Rn detector with the large PIN-photodiode is $2.32^{+0.34}_{-0.23}$~count/day while that of Rn detector with the small PIN-photodiode is $2.17^{+0.33}_{-0.22}$~count/day as summarized in Table~\ref{table_bg_setup}.

\begin{figure*}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{bg18_rate.eps}
\end{minipage}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{bg28_rate.eps}
\end{minipage}
\caption{The measured background rate of the Rn detectors as a function of time. The right (left) panel shows the result from the Rn detector with $18\times18~\mathrm{mm}$ ~($28\times28~\mathrm{mm}$) PIN-photodiode. The horizontal dashed line~(gray band) shows the average~(uncertainty) of background rate for each detector.}
\label{fig_bg_rate2}
\end{figure*}

\begin{table*}
\centering
\caption{Summary of the background rates of the calibration setup. We took $30$~days data for this measurement. The count rate shows the actual count rate including the intrinsic background from the two Rn detectors while the Rn concentration is obtained after subtracting the intrinsic backgrounds described in Sect.~\ref{sec_bg}.
\label{table_bg_setup}
\begin{tabular}{ccc
\hline
 &  Rn detector  & Rn detector  \\
 & with $28\times28~\mathrm{mm}$ & with $18\times18~\mathrm{mm}$ \\
\hline
Count rate~[count/day]& $2.32^{+0.34}_{-0.23}$ & $2.17^{+0.33}_{-0.22}$ \\
(Including intrinsic background & & \\
from two Rn detectors) & & \\ \hline
Rn concentration~[$\mathrm{mBq/m^{3}}$] & $0.40^{+0.19}_{-0.13}$ & $0.38^{+0.17}_{-0.13}$ \\
(Intrinsic backgrounds subtracted) & & \\
\hline
\end{tabular}
\end{table*

The background rates measured by two Rn detectors include their intrinsic backgrounds. Thus, the background from the calibration setup should be obtained by subtracted their intrinsic backgrounds from the measured background rates. For the Rn detector with $28\times28~\mathrm{mm}$ PIN-photodiode, the intrinsic background is the measurement result described in Sect.~\ref{sec_bg}. On the other hand, we have not performed the background measurement for the Rn detector with $18\times18~\mathrm{mm}$ PIN-photodiode in this study. To compensate for this deficiency of background data, we used the past result from Ref.~\cite{Nakano:2017rsy} to subtract the background rate. Finally, the measured Rn concentrations are $0.40^{+0.19}_{-0.13}~\mathrm{mBq/m^{3}}$ by $28\times28~\mathrm{mm}$ and $0.38^{+0.17}_{-0.13}~\mathrm{mBq/m^{3}}$ by $18\times18~\mathrm{mm}$ after subtracting the intrinsic backgrounds of two Rn detectors.

Using the correction factor~$\varepsilon$ defined in Eq.~(\ref{eq_3}), the daily Rn emanation rate is also estimated to be $0.04^{+0.03}_{-0.02}~\mathrm{\mu Bq/setup/day}$. Although some excesses over the intrinsic background are observed, the Rn emanation form the calibration setup is quit low. Moreover, there are no air leaks in the calibration setup since the environmental Rn concentration in the experimental site is about $40$--$80~\mathrm{Bq/m^{3}}$, which has been continuously monitored by the $1$~L Rn detector~\cite{Pronost:2018ghn}. Therefore, the background from the calibration setup does not affect the calibration results described in Sect.~\ref{sec_setup} because of its small contamination.

\section{Background measurement} \label{sec_bg}

\subsection{Background rate with purified air}
The intrinsic background rate was evaluated by closing the detector. This measurement evaluates the intrinsic background due to the emanation from the inner surface of the Rn detector, including the PIN-photodiode, because some amount of Rn emanating from other components of the calibration setup increases the background rate. To eliminate such uncertainties, we measured the background rate by closing the Rn detectors. In addition, the background from the calibration setup was also evaluated by circulating air. The result is described in Appendix~\ref{sec_bg_setup}.

Figure~\ref{fig_bg_spect} shows the measured daily background spectrum of the $80$~L Rn detector with $28\times28~\mathrm{mm^{2}}$ PIN-photodiode\footnote{The signal region is slightly different from that used in Figure~\ref{fig_range} since we replaced the electronics after the calibration described in Sect.~\ref{sec_setup}.} and the average background count rate\footnote{Once the Rn detector is closed, the residual Rn in the purified air starts to decay, and Rn also starts to emanate from the detector's inner surface. They contribute as backgrounds. To consider these two physics processes, it is proper to fit the background rate with the following function:
\begin{equation}
  Ae^{-\lambda t} + B(1-e^{-\lambda t}), \notag
\end{equation}
\noindent
where $t$ is the elapsed time since the start of the measurement, $\lambda$ is the decay constant for $\mathrm{^{222}Rn}$, and $A$ and $B$ are parameters to be fitted. However, the background rate shown in Figure~\ref{fig_bg_spect} is quite low. Hence, we calculated the average count rate as the background for this study.}. The count rate is $0.57^{+0.20}_{-0.10}$~count/day, which corresponds to a Rn concentration of $0.24^{+0.09}_{-0.05}~\mathrm{mBq/m^{3}}$. Table~\ref{table_bg} summarizes the background rates evaluated in this study, as well as those measured in previous studies. The intrinsic background of the $80$~L Rn detector with $28\times28~\mathrm{mm}$ is consistent with that of the $18\times18~\mathrm{mm}$ evaluated in previous studies~\cite{Hosokawa:2015koa, Nakano:2017rsy}. Moreover, their background rates were unchanged, although we installed a larger PIN-photodiode.

\begin{figure*}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{closed_bg_spect.eps}
\end{minipage}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{closed_bg_rate.eps}
\end{minipage}
\caption{(Left) The measured background spectrum of the $80$~L Rn detector with $28\times28~\mathrm{mm}$ of PIN-photodiode. The vertical axis shows the count/channel/day, and the bottom~(top) horizontal axis shows the ADC channel~(measured energy in units of MeV). The two dashed lines show the upper and lower bounds of the $\mathrm{^{214}Po}$ signal region. (Right) The $\mathrm{^{214}Po}$ count rate of the 80~L Rn detector with $28\times28~\mathrm{mm}$ of the PIN-photodiode is a function of time. The blue square points show the daily count rate, and the dashed horizontal line~(gray band) shows the average~(uncertainty) of the count rate.}
\label{fig_bg_spect}
\end{figure*}

\begin{table*}
\centering
\caption{Summary of the intrinsic background rates of the $80$~L Rn detector equipped with the $28\times28~\mathrm{mm}$~($18\times18~\mathrm{mm}$) PIN-photodiode. These measurements were conducted by filling the purified air, and the absolute humidity was controlled to be less than $10^{-3}~\mathrm{g/m^{3}}$.}
\label{table_bg}
\begin{tabular}{cccc}
\hline
 & Livetime & Rn detector & Rn detector  \\
 & [day] &  with $28\times28~\mathrm{mm}$ & with $18\times18~\mathrm{mm}$ \\
 \hline
 This study & $28$ & $0.57^{+0.20}_{-0.10}$~count/day & --  \\
 ($80$~L Rn detector) & & $0.24^{+0.09}_{-0.05}~\mathrm{mBq/m^{3}}$ & --   \\
\hline \hline
 Past study 1~\cite{TAKEUCHI1999334} & Approximately  &  -- &  $2.4\pm1.3$~count/day \\
 ($70$~L Rn detector) & 2~months & -- & -- \\
\hline
 Past study 2~\cite{Hosokawa:2015koa} & $132.5$ &  -- &  $0.81\pm0.08$~count/day \\
 ($80$~L Rn detector) & & -- & $0.37\pm0.05~\mathrm{mBq/m^{3}}$ \\
\hline
 Past study 3~\cite{Nakano:2017rsy} & $156$ &  -- &  $0.78\pm0.07$~count/day \\
 ($80$~L Rn detector) & & -- & $0.33\pm0.07~\mathrm{mBq/m^{3}}$ \\
\hline
\end{tabular}
\end{table*}

Relatively large signals from $\mathrm{^{210}Po}$ are observed, as shown in Figure~\ref{fig_bg_spect}. These signals come from the $\mathrm{^{222}Rn}$ daughters remain on the PIN-photodiode once they have been collected because the intrinsic background measurement was conducted after the calibrations described in Sect.~\ref{sec_setup}.

The daily Rn emanation rate from the $80$~L Rn detector can be estimated by considering the correction of the decay duration since its generation. The correction factor, defined as $\varepsilon$, is expressed as
\begin{align}
\varepsilon & = \frac{\mathrm{Number~of~Rn~after~emanation}}{\mathrm{Number~of~Rn~under~radioactive~equilibrium}} \notag \\
 & = \frac{N\left\{1-\left(\frac{1}{2}\right)^{1.0/3.82}\right\} }{N\left\{1-\left(\frac{1}{2}\right)^{\infty/3.82}\right\} } = 0.1659. \label{eq_3}
\end{align}

\noindent Using this correction factor, the daily emanation rate from the $80$~L Rn detector is calculated as $(3.18^{+1.19}_{-0.66})~\mathrm{\mu Bq/detector/day}$.

\section{Calibration with purified air} \label{sec_setup}
\subsection{Calibration setup}
We constructed a calibration setup at Lab-A in the experimental area at Kamioka Observatory, the Institute for Cosmic Ray Research~(ICRR), and the University of Tokyo. Figure~\ref{fig_setup} shows the schematic view of the calibration setup and the image of the actual setup.

\begin{figure*}
\centering\includegraphics[width=\linewidth]{setup.eps}
\caption{Schematic view of the calibration setup~(left) and image of the actual setup~(right). Arrows in the left schematic view show the direction of air circulation.}
\label{fig_setup}
\end{figure*}

The system consists of two Rn detectors, a temperature sensor~(VISALA DMT 340), an air circulation pump~(IWAKI BA-330SN), a cold trap~(Taisho TC0147), and pressure gauges~(Naganokeiki Co. Ltd. ZT67), and a mass flow controller~(HORIBA SEC-Z512MGX). We installed a large PIN-photodiode to one of the two Rn detectors and a small PIN-photodiode to the other. To prevent dust from entering the Rn detector, a $3$~nm filter~(Pall Corp. GLFF4281VMM4) was installed in front of each Rn detector. To supply a high voltage to the Rn detectors in parallel, a high-voltage divider~(HAYASHI-REPIC Co., Ltd., RPH-033) was used. To reduce the electrical noise, a noise cut trans~(DENKENSEIKI Research Institute Co., Ltd., NCT-I3) was connected.

The temperature sensor monitors the temperature and dew point of the circulating gas. The pressure gauge was also equipped with a Rn detector to monitor the inner pressure. To control the humidity in the gas, we used a cold trap, which is a 1/2-inch U-shaped pipe filled with $14.0$~g of $\phi~80~\mu\mathrm{m}$ copper wool~(Nippon Steel Wool Co. Ltd.). This U-shaped pipe was placed in a heat-insulated vacuum container and connected to a refrigerator.

Other system components are made of electro-polished stainless steel~(NISSHO Astec Co. Ltd. MGS-EP SUS316L), and all joints are connected by VCR$^{\textcircled{\footnotesize{R}}}$ gaskets to minimise Rn emanation and possible air leaks, both of which may affect the background measurements. The air leak rate of the system was measured to be less than $10^{-10}~\mathrm{Pa\cdot m^{3}/s}$ using a helium leak detector~(ULVAC Equipment Sales Inc. HELIOT 712D2).

To inject Rn gas into the calibration system, we used an Ra source~(PYLON RNC) whose activity is ($199\pm4\%$)~Bq. This calibration source was connected to the system with an adapter of a steel quick connect stem produced by Swagelok, as shown in Figure~\ref{fig_setup}.

\subsection{Calibration procedure}
Calibrations were performed by the following procedures. Before the calibration, we evacuated the entire calibration setup down to $1.0\times10^{-4}$~Pa by a vacuum pump~(Hakuto Co. Ltd. HiCube 80 Eco). Then, we filled the entire calibration system with commercially available G1-grade purified air, whose impurities were less than $0.1$~ppm, until the inner pressure reached the atmospheric pressure. We attached the Ra source to the calibration system and started to circulate the air to maintain a constant Rn concentration in the calibration system. After reaching radioactive equilibrium in the calibration system, we measured the count rate of the $\mathrm{^{214}Po}$ signal from the Rn detector.

The Rn concentration under radioactive equilibrium was stable at approximately $1243.8\pm49.7~\mathrm{Bq/m^{3}}$, where the total volume of the calibration setup is $0.160~\mathrm{m^{3}}$. Note that the total volume, except for the two Rn detectors, is only $0.6\times10^{-3}~\mathrm{m^{3}}$, and this volume can be ignored in the calculation because the uncertainty of the Rn detector volume is approximately $2$--$3\times10^{-3}~\mathrm{m^{3}}$ level, as mentioned in the later section~(Sect.~\ref{sec_sys}). Finally, we evaluated the calibration factors~($C_{F}$) by comparing the count rate and constant Rn concentration.

\subsection{High-voltage dependence}
Because the electrical field between the PIN-photodiode and the stainless steel vessel is formed by the supplied high voltage, the collection efficiency is highly dependent on the supplied high-voltage value. To confirm this dependence, we evaluated the calibration factor by changing the supplied high voltage. Figure~\ref{fig_spect_shift} shows examples of the daily spectrum measured by the Rn detector equipped with the small~(large) PIN-photodiode during the calibration.

\begin{figure*}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{spect18.eps}
\end{minipage}
\begin{minipage}{0.5\hsize}
\centering\includegraphics[width=1.0\linewidth]{spect28.eps}
\end{minipage}
\caption{Typical daily spectrum measured by the Rn detector equipped with  $18\times18~\mathrm{mm}$~($28\times28~\mathrm{mm}$) PIN-photodiode in left (right). As an example, the spectrum shown here was measured under radioactive equilibrium by supplying high voltage at $-2.0$~kV~(red in left and blue in right) and $-0.2$~kV~(yellow in left and green in right).}
\label{fig_spect_shift}
\end{figure*}

There are two large peaks in the spectrum originating from $\mathrm{^{218}Po~(6.00~MeV)}$ and $\mathrm{^{214}Po~(7.69~MeV)}$. The $\mathrm{^{214}Po}$ signal is used to calculate the Rn concentration because there are no other $\alpha$~sources near its peak energy\footnote{In the case of the $\mathrm{^{218}Po}$ signal, several background sources from the $\mathrm{^{232}Th}$ series exist. In particular, $\alpha$ particles are emitted from $\mathrm{^{208}Tl}~(5.00~\mathrm{MeV})$,  $\mathrm{^{228}Th}~(5.52~\mathrm{MeV})$, $\mathrm{^{224}Ra}~(5.79~\mathrm{MeV})$, $\mathrm{^{212}Bi}~(6.21~\mathrm{MeV})$, and $\mathrm{^{220}Rn}~(6.41~\mathrm{MeV})$. Their signals potentially overlapped with the peak energy of $\mathrm{^{218}Po}$. Thus, $\mathrm{^{218}Po}$ signals are not used to calculate the Rn concentration in the gases.}. Even though the PIN-photodiode has a high spectral resolution, background signals from $\mathrm{^{216}Po~(6.78~MeV)}$ and $\mathrm{^{212}Po~(8.95~MeV)}$ potentially overlap with the signals from $\mathrm{^{214}Po}$. Hence, the signal window must be properly determined to evaluate the calibration factor. For this reason, we set the signal window to count $\mathrm{^{214}Po}$ signals based on the criteria in a previous publication~\cite{Nakano:2017rsy}. We set the lower bound at five channels from the peak of $\mathrm{^{216}Po}$ signals not only to prevent the entry of the $\mathrm{^{216}Po}$ signals but also to cover the tail of the $\mathrm{^{214}Po}$ signals. On the other hand, we set the upper bound at five channels from the peak of $\mathrm{^{214}Po}$ signals to prevent the entry of $\mathrm{^{212}Po}$ signals. Note that $\mathrm{^{216}Po}$ and $\mathrm{^{212}Po}$ originate from the decay of $\mathrm{^{220}Rn}$, and their observed rates are quite low compared with that of $\mathrm{^{214}Po}$ signals during the calibration.

The pulse height of the output signal shifts when the supplied high-voltage changes because the electronics simultaneously supply the reverse bias minus the high voltage to the PIN-photodiode. Figure~\ref{fig_spect_shift} also shows the daily spectrum by changing the supplied high voltage from $-0.2$~kV to $-2.0$~kV. As expected, the peak positions of the observed $\mathrm{^{214}Po}$~(and $\mathrm{^{218}Po}$) signals shift depending on the supplied high voltage. Furthermore, the total number of daily observed signal rates decreases when the supplied high voltage is lowered.

\begin{figure}
\centering\includegraphics[width=\linewidth]{spect_range.eps}
\caption{The region of the ADC channel determined from the peak position of $\mathrm{^{214}Po}$ signals observed by the PIN-photodiode as a function of the supplied high voltage~[kV]. The red circle~(blue square) points show the ADC channel of the signal peak position for the PIN-photodiode of $18\times18~\mathrm{mm}$~($28\times28~\mathrm{mm}$). The red~(blue) bands show the windows for selecting $\mathrm{^{214}Po}$ signals based on the criteria used in a previous study~\cite{Nakano:2017rsy}.}
\label{fig_range}
\end{figure}

Figure~\ref{fig_range} shows the selected signal windows for the calibration as a function of the supplied high voltage. The peak position of $\mathrm{^{214}Po}$ signals relatively shifts significantly below $-1.0$~kV while the peak position does not shift considerably above $-1.0$~kV.

By using these signal windows, we obtained the calibration factors for each high-voltage value. Figure~\ref{fig_hv_dep} shows the calibration factors as functions of the supplied high voltage. To quantitatively evaluate the difference between the two PIN-photodiodes, the relative ratio between their calibration factors, which is defined as
\begin{equation}
    R_{C_{F}} = \frac{C_{F,\,\mathrm{large}}-C_{F,\,\mathrm{small}}}{C_{F,\,\mathrm{small}}} \times 100.0~[\%], \label{eq_2}
\end{equation}

\noindent
is also calculated for each high-voltage setting. Based on the measurements, we found a $15$--$30\%$-level of improvement below $-1.0$~kV, while only a few~$\%$ level of change occurred above this level. This behaviour is similar to that observed in the region where the electrical field is formed, as shown in Figure~\ref{fig_active}.

We chose $-2.0$~kV for the calibration of absolute humidity as well as the background rate measurement in the later sections because the divided high-voltage supply to the PIN-photodiode becomes larger than its tolerance limit. Note that the calibration factor is almost saturated above $-2.0$~kV based on the simulation shown in Figure~\ref{fig_active}, and this setting does not significantly affect the calibration result.

\begin{figure}
\centering\includegraphics[width=\linewidth]{cf_hv.eps}
\caption{The supplied high-voltage dependence of the calibration factor. In the top panel, the blue square~(red circle) points show the calibration factor of the Rn detector equipped with a $28\times28~\mathrm{mm}$~($18\times18~\mathrm{mm}$) PIN-photodiode determined by changing the supplied high voltage. The calibrations were performed by setting the air circulation rate to $0.6$L/min. To ease visualisation, the vertical error bars were scaled by $2$, and the red circles horizontally shifted by $-0.25$~kV from the original position. The absolute humidity during this calibration was controlled at $0.004~\mathrm{g/m^{3}}$. In the bottom panel, the black triangle points show their relative ratio~($R_{C_{F}}$) defined in Eq.~(\ref{eq_2}).}
\label{fig_hv_dep}
\end{figure}

\subsection{Absolute humidity dependence} \label{sec_ah}
The collection efficiency also depends on the humidity of the gas because the charged nuclei are captured by water. Their neutralisation results in a loss of collection efficiency~\cite{osti_5305900, howard1991}. To test this dependence, we calibrated $C_{F}$ by changing the temperature of the cold trap from $-80\mathrm{^{\circ}C}$ to room temperature~(typically $+20\mathrm{^{\circ}C}$ in the Kamioka mine).

Figure~\ref{fig_ah_dep} shows the calibration factors as a function of the absolute humidity of filling purified air and the relative ratio between them~($R_{C_{F}}$). By comparing the calibration factor of the small PIN-photodiode with the large one, we found a level of improvement of $10$--$20\%$ above $1.0~\mathrm{g/m^{3}}$ while a $(3.8\pm2.4)\%$-improvement was observed below this level by replacing the PIN-photodiode.

\begin{figure}
\centering\includegraphics[width=\linewidth]{cf_ah.eps}
\caption{The absolute humidity dependence of the calibration factor. In the top panel, the blue square~(red circle) points show the calibration factor determined by the Rn detector equipped with a $28\times28~\mathrm{mm}$~($18\times18~\mathrm{mm}$) PIN-photodiode. The calibrations were performed by setting a high voltage of $-2.0$~kV and an air circulation rate of $0.6$~L/min. To ease visualisation, the vertical error bars are scaled by $2$, and the red circles horizontally shift by $-10\%$ from the original position. In the bottom panel, the black triangle points show the relative ratio~($R_{C_{F}}$) defined in Eq.~(\ref{eq_2}).}
\label{fig_ah_dep}
\end{figure}

The behaviour of this dependence is well described by an empirical function, $C_{F}=A-B\sqrt{A_{H}}$, where $A$ and $B$ are parameters to be fitted, and $A_{H}$ is the absolute humidity of the filling gas~\cite{osti_6605848}. Figure~\ref{fig_ah_band} shows a comparison of the empirical function obtained in this study as well as those of previous studies~\cite{TAKEUCHI1999334, Nakano:2017rsy}, where the fitting parameters are summarized in Appendix~\ref{sec_app1}.

\begin{figure}
\centering\includegraphics[width=\linewidth]{comp_cf_band.eps}
\caption{Comparison of the absolute humidity dependence of the calibration factor among the Rn detectors when the purified air is filled in. The blue and red lines~(with band) show the empirical functions obtained in this study. The green, black, and pink lines show the results of previous studies~\cite{TAKEUCHI1999334, Nakano:2017rsy}.}
\label{fig_ah_band}
\end{figure}

In the case of the $70$~L Rn detector, the acrylic plate and Viton O-rings were used as the top flanges of the detector. Because of this design, the detector is not airtight, and therefore its sensitivity is worse in humid conditions~\cite{TAKEUCHI1999334}. After replacing the top acrylic flange with knife-edge flanges with metal gaskets, the calibration factors in humid conditions were improved~\cite{Hosokawa:2015koa, Nakano:2017rsy}.

By comparing the shapes of the empirical functions, we found that the $80$~L Rn detector with the PIN-photodiode of $28\times28~\mathrm{mm}$ has the best sensitivity among them. This improvement under dry conditions below $0.1~\mathrm{g/m^{3}}$ results from replacing the PIN-photodiode.

\subsection{Systematic uncertainty in the calibration} \label{sec_sys}
The systematic uncertainties are estimated according to the technical specifications of the measurement devices. The assigned systematic uncertainties are summarized in Table~\ref{table_sys}.

\begin{table}
\caption{Summary of the systematic uncertainties in the calibration.}
\label{table_sys}
\centering
\begin{tabular}{cc}
\hline
Item & Assigned \\
 & uncertainty \\
\hline
Ra source activity& $\pm4\%$\\
Temperature sensor& $\pm2\%$\\
Air flow rate& $\pm2\%$\\
Volume of the whole calibration setup& $\pm3\%$\\
\hline
\end{tabular}
\end{table}

Based on the manufacturer's documentation, the uncertainty of the activity of the Ra source is $\pm4\%$. The temperature sensor has an uncertainty of $\pm2\%$, and the mass flow controller has an uncertainty of $\pm2\%$. Because the volume of the whole calibration system is not measured precisely, we assigned the systematic error of the total volume as $\pm3\%$, where the dominant uncertainty comes from the accuracy of the Rn detector's structure. These systematic uncertainties are included in the calibration results described in this section.

\section{Summary and conclusion} \label{sec_conclusion}
To improve the detection sensitivity for measuring Rn in gases, we developed a new Rn detector by replacing the PIN photodiode, HAMAMATSU-S3584-09. To evaluate the detector performance, we constructed a calibration setup by serially connecting two Rn detectors equipped with both types of PIN photodiodes.

First, we tested the high-voltage dependence of the calibration factor. Based on the calibrations, we found a $15$--$30\%$-level of improvement of the calibration factor below $-0.3$~kV, while a few~$\%$ improvements occurred above this level. We also tested the absolute humidity dependence of the calibration factor. We found an improvement of $(3.8\pm2.4)\%$ in the detection efficiencies below $1.0~\mathrm{g/m^{3}}$, while a $10$--$20\%$-improvement occurred above this level. We also measured the intrinsic background of the detector, $0.24^{+0.09}_{-0.05}~\mathrm{mBq/m^{3}}$. This background is consistent with the Rn detector of the small PIN-photodiode. We successfully demonstrated that the sensitivity of the detector was improved by replacing the larger PIN-photodiode.

The limitation of the collection efficiency in the current design comes from the difficulty in forming the electrical field between the PIN-photodiode and the stainless steel structure because of its small volume. Further sensitivity is expected by enlarging the Rn detector's volume.

\section{Rn detector} \label{sec_detector}
\subsection{Detector structure and electronics}
The Rn detector consists of a stainless steel vessel, a PIN-photodiode, a ceramic feed through, and an electrical circuit. Figure~\ref{fig_det} shows the structure of the $80$~L Rn detector and the electrical field simulated by COSMOL Multiphysics$^{\textcircled{\footnotesize{R}}}$\footnote{\url{https://www.comsol.jp/comsol-multiphysics}}. The PIN-photodiode was connected to the inner end of the feedthrough, which was mounted on the top flange of the vessel. The high-voltage divider and readout electronics are also connected to the external end of the feedthrough. The details of the Rn detector can be found in Ref.~\cite{Hosokawa:2015koa, Nakano:2017rsy}.

\begin{figure*}
\centering\includegraphics[width=0.9\linewidth]{det_sim_revised2.eps}
\caption{Left: The conceptual diagram of the Rn detector and the electrical field simulated by COMSOL Multiphysics$^{\textcircled{\footnotesize{R}}}$. In this simulation, $28\times28~\mathrm{mm}$ of the PIN-photodiode is installed at the top of the detector, and $-2.0$~kV is supplied to form the electrical field. The red solid lines show the equipotential lines~($-2000$~V, $-1000$~V, $-500$~V, $-100$~V, $-50$~V, $-10$~V, $-5$~V, and $-1.0$~V, respectively). Right: Image of the Rn detector and the PIN-photodiode.}
\label{fig_det}
\end{figure*}

The electronics consisted of a high-voltage supplier and an amplifier circuit. The reverse bias voltage was simultaneously supplied to the PIN-photodiode. In this study, we used a Raspberry Pi-based system originally developed for monitoring the Rn concentration in the underground experimental area in the Kamioka mine with a $1$~L Rn detector~\cite{Pronost:2018ghn}. The data acquisition system recorded the signal from the Rn detector and the environmental parameters, such as the temperature, dew point, inner pressure of the detector, and flow rate of gas circulation in the calibration system.

\subsection{Larger size PIN-photodiode and formed electrical field}
In previous studies, we used various sizes of PIN-photodiodes produced by Hamamatsu Photonics K. K., such as S3590-06~($10\times10~\mathrm{mm}$)~\cite{1994125}, S3590-09~($10\times10~\mathrm{mm}$)~\cite{Pronost:2018ghn}, S3204-06~($16\times16~\mathrm{mm}$)~\cite{1997710, TAKEUCHI1999334, Mitsuda:2003eu}, and S3204-09~($18\times18~\mathrm{mm}$)~\cite{Hosokawa:2015koa, Nakano:2017rsy}.  In this study, S3584-09 with a sensitive area of $28\times28~\mathrm{mm}$ was used. Its specifications can be found on the company's website\footnote{\url{https://www.hamamatsu.com/resources/pdf/ssd/s3204-08_etc_kpin1051e.pdf}}.
Figure~\ref{fig_pin} shows an example of the PIN-photodiodes used in this study. The new PIN-photodiode has approximately $2.5$ times larger surface area than that of previously used PIN-photodiode~(S3204-09, $18\times18~\mathrm{mm}$). In this article, we refer the PIN-photodiode with $18\times18~\mathrm{mm}$ surface area to the small PIN-photodiode while that with $28\times28~\mathrm{mm}$ to the large PIN-photodiode.

\begin{figure}
\centering\includegraphics[width=0.8\linewidth]{pin_comp.eps}
\caption{A picture of PIN-photodiodes used in this study. The left is HAMAMATSU-S3204-09~($18\times18~\mathrm{mm}$), while the right is HAMAMATSU-S3584-09~($28\times28~\mathrm{mm}$).}
\label{fig_pin}
\end{figure}

As briefly explained in Sect.~ \ref{sec_intro}, the volume where the electrical field is formed increases when the PIN-photodiode with a larger size is installed. Before the calibration described in Sect.~\ref{sec_setup}, we simulated the volume where the electrical field is formed with COSMOL Multiphysics$^{\textcircled{\footnotesize{R}}}$. In this simulation, we assumed a current detector structure with a PIN-photodiode at the top. To quantitatively evaluate the difference, we define the proportion of the volume where the electrical field is formed above $-1.0$~V to the total inner volume of the $80$~L Rn detector as $P_{\mathrm{small}}~(P_{\mathrm{large}})$ when the $18\times18~\mathrm{mm}$~($28\times28~\mathrm{mm}$) PIN-photodiode is installed on the Rn detector.

Figure~\ref{fig_active} shows $P_{\mathrm{small}}$ and $P_{\mathrm{large}}$ as functions of the supplied high voltage and their relative ratio, which is defined as
\begin{equation}
R = \frac{P_{\mathrm{large}}-P_{\mathrm{small}}}{P_{\mathrm{small}}} \times 100.0~[\%]. \label{eq_1}
\end{equation}
The volume where the electrical field is formed proportionally increases to the supplied high voltage. The relative ratio~($R$) rapidly increases until the supplied high voltage reaches $-2.0$~kV, while it is almost saturated above this level. For example, $P_{\mathrm{small}}$ is $85.63\%$, whereas  $P_{\mathrm{large}}$ is $88.71\%$ at $-2.0$~kV. The difference between them is $+3.59\%$ at $-2.0$~kV, which demonstrates that the electrical field formed inside the detector is enlarged by replacing the PIN-photodiodes. Therefore, the large PIN-photodiode allows the Rn detector to improve its sensitivity. Below $-2.0$~kV, the difference between $P_{\mathrm{samll}}$ and $P_{\mathrm{large}}$ is significantly large, which implies that the small PIN-photodiode is screened by the surrounding stainless structure.

\begin{figure}
\centering\includegraphics[width=\linewidth]{hv_ratio_revised.eps}
\caption{Top: The proportion~($P_{\mathrm{small}}$ or $P_{\mathrm{large}}$) of the volume where the electrical field is formed above $-1.0$~V to the total inner volume of the $80$~L Rn detector. The vertical(horizontal) axis shows the ratio in units of percent~(supplied high voltage in units of kV). The blue square~(red circle) points show the simulated result when the PIN-photodiode of $28\times28~\mathrm{mm}$~($18\times18~\mathrm{mm}$) is installed on the $80$~L Rn detector. For ease of visualisation, the red circle shifts horizontally by $-0.25$~kV from the original position. Bottom: The black triangle shows the relative ratio~($R$) defined in Eq.~(\ref{eq_1}).}
\label{fig_active}
\end{figure}

\subsection{Calibration factor}
Regarding the collection efficiency of the Rn detector, we defined a calibration factor~($C_{F}$) as a parameter to convert the count rate of $\mathrm{^{214}Po}$ in units of count/day into the Rn concentration in units
of $\mathrm{mBq/m^{3}}$. In the calibration described in the next section, we evaluated $C_{F}$ by measuring the count rate of $\mathrm{^{214}Po}$ signals from a particular Rn detector with a constant Rn concentration under radioactive equilibrium in a calibration setup.

\section{Application of the Rn detector with a large PIN-photodiode} \label{sec_future}
\subsection{Rn measurement of purified water}

As prospects, upgrading by replacing the larger PIN-photodiode is helpful to precisely monitor its Rn concentration in a humid gas, as described in Sect.~\ref{sec_ah}. In the case of Super-Kamiokande, the Rn concentration in the purified water was measured by the Rn detector combined with a membrane degasifier module~\cite{TAKEUCHI1999334, Mitsuda:2003eu}. However, the sensitivity of this method is limited by the inefficiency of the Rn detector in a humid gas after extracting Rn gas from the purified water. Another approach using both an extraction column and a cooled charcoal trap is adapted at the expense of real-time measurement instead of the membrane degasifier module~\cite{Nakano:2019bnr}. By replacing the large PIN-photodiode and the newly developed membrane degasifier module with a low intrinsic background~\cite{takeda:nu2020}, the sensitivity of the Rn detector combined with the membrane degasifier module is expected to be improved.

\subsection{Material screening by measuring the Rn emanation rate}

In rare event search experiments, the radioactive impurities in the detector's material are severe background sources. The emanation rate is fixed once the detector is constructed. For this reason, careful material screening programs with several screening devices have been conducted before the construction of each detector ~\cite{Leonard:2007uv, Aprile:2011ru, Wang:2016eud, Aprile:2017ilq, Cebrian:2017jzb}.

The Rn detector is used to measure the Rn emanation rate from materials by combining it with a sampling chamber in which the material is put in~\cite{Sekiya:2015rla,Takeuchi:2017pij}. In such a measurement system, water is also emitted from the surface of the material, resulting in the neutralisation of Rn daughters. Hence, the control of absolute humidity is the key to precisely evaluating the material's Rn emanation rate. In this study, we found a higher collection efficiency in humid gases by replacing a small PIN-photodiode with a larger one. Therefore, this upgrade also helps to apply the Rn detector to such screening programs to evaluate the impurities in the detector's materials.

\section{Introduction} \label{sec_intro}

Noble gas radon~($\mathrm{^{222}Rn}$, hereinafter Rn) is one of the major background sources in experiments that search for rare~(interaction/decay) events, such as solar neutrinos, neutrino-less double beta decay, and dark matter. Once a detector for a particular experiment is constructed, Rn gas is continuously generated from the decay chain of the $\mathrm{^{238}U}$ series in the detector material. Owing to its long half-life~($3.82$~days), its diffusion into the detector's fiducial volume results in a mimicking signal in that energy range. Therefore, material screening before the detector construction, monitoring Rn in target materials used in such experiments, and Rn removal during the detector operation by their purification processes are essential for improving their sensitivities.

For many years, several technologies for Rn measurements have been developed for these purposes. In particular, a Rn detector equipped with a PIN-photodiode is widely used to measure the Rn concentration in gases~\cite{4329405, HOWARD1990589, 1994125, 1996741, 1997710, TAKEUCHI1999334, CHOI2001177, KIKO2001272, Mitsuda:2003eu, GUTIERREZ2004583, MARTINMARTIN20061287, Mamedov:2011zz, 2013InJPh..87..471A, Hosokawa:2015koa, Nakano:2017rsy, Pronost:2018ghn, Elisio:2019nuz}. The principal technique of such Rn detectors is electrostatic collection~\cite{Kotrappa1981ElectretaNT}. Because the daughter nuclei of $\mathrm{^{222}Rn}$ tend to be positively charged~\cite{PMID:2663781}, they are easily collected by forming an electrical field inside the detector. The PIN-photodiode, where charged nuclei are collected, uniquely identifies the energy of the $\alpha$ particle emitted via the daughter's subsequent decays because of its high spectral resolution.

For the past $30$~years, we have developed this type of Rn detector. The current design of this detector, whose inner volume is approximately $80$~litter~(hereinafter~L), has been used to measure Rn concentrations less than a few $\mathrm{mBq/m^{3}}$ in several gases~\cite{Hosokawa:2015koa, Nakano:2017rsy}. Recently, we applied this detector to measure the Rn concentration in purified water in the Super-Kamiokande detector by combining it with a Rn extraction column and cooled charcoal~\cite{Nakano:2019bnr} and to evaluate several absorbents for removing Rn from purified gases~\cite{Ogawa:2019ccp, Nakano:2020tej}.

For future large-scale experiments to search for rate events, further sensitivity is required to monitor the Rn concentration below the $\mu\mathrm{Bq/m^{3}}$ level in purified gases~\cite{Abe:2011ts, Aalseth:2017fik, Aprile:2020vtw, Agostini:2020adk}. There are four strategies for improving the detection sensitivity of such Rn detectors. The first is to enlarge the detector volume, the second is to reduce the intrinsic background, the third is to raise the supplying high voltage to form a large electric field inside the detector, and the fourth is to enlarge the surface area of the PIN-photodiode. The details of these strategies are as follows.

\begin{description}
\item[1.~Detector volume] \mbox{}\\
The enlargement of the detector volume is an easy way to improve the sensitivity. For instance, to measure Rn concentrations less than a few $\mathrm{mBq/m^{3}}$ in gases, several large Rn detectors have been developed in the Borexino collaboration~($480$~L~\cite{Simgen:2003et}), Super-Kamiokande collaboration~($700$~L~\cite{Mitsuda:2003eu}), and JUNO collaboration ($279$~L~\cite{Yu:2020vbz}). On the other hand, this strategy has two disadvantages: it is cost-intensive and increases the intrinsic background emanating from the inner surface of the detector structure.
\item[2.~Intrinsic background] \mbox{}\\
The background of such a detector comes from Rn emanating from the inner surface as well as air leaks. Electropolishing the inner surface efficiently reduces Rn emanating from detector materials~\cite{TAKEUCHI1999334}. Furthermore, the latest $80$~L detector used knife-edge flanges with metal gaskets~(CF flanges) to prevent the entry of external air~\cite{Hosokawa:2015koa, Nakano:2017rsy}. This technique enables the detector to be airtight. Based on these techniques, the background rate of the latest $80$~L Rn detector was evaluated as $0.74\pm0.07~\mathrm{(Statistical~uncertainty~only)}$~count/day when purified air was filled with~\cite{Nakano:2017rsy}. Owing to its low background rate, the sensitivity was less than $0.5~\mathrm{mBq/m^{3}}$ level by a single-day measurement.

Another possibility is copper electroplating because the radio-contaminants of copper materials are generally low~\cite{hoppe2008, Alduino:2016vjd, Abe:2017jzw, Bunker:2020sxw}. This technique may reduce the intrinsic background from the inner surface of the detector because the metal layer on its surface acts as an internal shield~\cite{Hoppe:2014nva, Balogh:2020nmo}. However, this approach, as well as the enlargement of the detector volume, is expensive.

\item[3.~High voltage] \mbox{}\\
High voltage directly affects the collection efficiency in principle because the charged nuclei are collected by the electrical field formed between the PIN-photodiode and a stainless steel structure. Therefore, increasing the supplied voltage is an easy way to improve the sensitivity of the detector. However, this approach depends on the shape of the detector and the location of the installed PIN-photodiode. In other words, the sensitivity of the detector is limited by the detector's size because of the difficulty in forming the electrical field efficiently inside the detector.

\item[4.~PIN-photodiode] \mbox{}\\
The surface area of the PIN-photodiode is sensitive to the collection efficiency because the charged nuclei are eventually collected in that area. Furthermore, the electrical field formed by the small size of the PIN-photodiode is screened by the surrounding stainless structure~\cite{2011JInst...6C2018F} as the electrical field formation is determined by an electrical dipole. Therefore, a larger surface area enables the Rn detector to more efficiently collect the charged nuclei and increase the volume where the electrical field is applied inside the detector.
\end{description}

Based on these discussions, we chose a PIN-photodiode whose surface area is $2.5$~times larger than that used previously as a first step. This approach enables the detector to improve its sensitivity at a low cost, while other approaches are cost-intensive.

This article consists of six sections, including an introduction. In Sect.~\ref{sec_detector}, we briefly introduce the developed Rn detector and describe its characteristics. In addition, we discuss the advantages of installing a large PIN-photodiode to the Rn detector. In Sect.~\ref{sec_setup}, we describe the setup of and results from the calibration using purified air. In particular, we report the high voltage and absolute humidity dependence of the collection efficiency and compare the results with those obtained in previous studies. In Sect.~\ref{sec_bg}, we evaluate the background of the Rn detector and compare it with past results. In Sect.~\ref{sec_future}, we discuss the applications of the Rn detector in precise measurement programs. Finally, we summarize and conclude the study in Sect.~\ref{sec_conclusion}.

\section*{Acknowledgment}

We gratefully acknowledge the cooperation of the Kamioka Mining and Smelting Company. We thank A.~Takeda from Kamioka Observatory, ICRR, the University of Tokyo, for preparing the additional electrical supply line at the experimental site. K. O. and Y. N. thank T. Onoue for constructing the support structures for the calibration setup. This work was partially supported by MEXT KAKENHI Grant Numbers 18H05536, 19J21344, 20K03998, and 21K13942. This work was partially supported by an inter-university research program at ICRR.

\bibliographystyle{ptephy}


\section{Introduction}
Learning quantum properties of an unknown physical system is essential for designing and testing devices based on the laws of quantum mechanics. Complete information about the physical state~\cite{QST,Gross2010} and process~\cite{Brien2004,Mohseni2008} of such a device can be obtained through tomography which requires considerable resources for implementation. On the other hand, some particular properties of a quantum system, like certification of quantum entanglement or incompatibility of measurements, can be learnt with less resources by constructing suitable witness operators~\cite{Guhne2009,Horodecki2009,Carmeli2019,Bae2020,Bavaresco2018}. The aforesaid methods still might not be optimal or available in some scenarios \cite{Reichardt2012,Vazirani2014,Miller2016}. Interestingly, however, some physical systems can be certified by employing comparatively far less resources through \textit{device-independent} tests where a device is treated simply as a black box~\cite{Scarani2012,Pironioetal2016, AcinandNavascues2017}. Then only from the input-output statistics termed as \textit{correlation} or \textit{behavior} of the box one can find the quantum state of the device. Such a certification is referred to as a self-test as it enables a user to verify the device without knowing details of its inner-workings~\cite{Mayers2004,McKague2011,YN2013,Bamps-Pironio2015,Wang-Wu-Scarani2016,Coladangeloetal,Coopmanset.al.2019,Jed2020,Ishizaka2020,Rabelo2012,DiCabello2021,Baccarietal2020,SupicBowles2020}.

Bell inequalities naturally fit into the device-independent paradigm for certification of quantum systems since their derivations, based on the assumptions of local-realism~\cite{EPR1935}, are independent from a quantum description of the physical state for the system or measurements applied to it~\cite{Bell1964,CHSH,Hardy1993,MJWHall2011,Scarani2019,Brunneretal2014}.
Therefore, various kind of Bell inequalities play a central role in the construction of self-testing protocols~\cite{Mayers2004,McKague2011,YN2013,Bamps-Pironio2015,Wang-Wu-Scarani2016,Coladangeloetal,Coopmanset.al.2019,Jed2020,Ishizaka2020,Rabelo2012,DiCabello2021,Baccarietal2020,SupicBowles2020}. In a Bell test only some special type of input-output statistics can self-test a quantum state; the one which can be realized (up to local isometry) by a unique quantum state and measurements~\cite{SupicBowles2020}.
Further, within the (convex) set of all quantum correlations in a specified Bell scenario, behaviors leading to self-testing protocols must necessarily be  extremal points of the quantum set~\cite{Rabelo2012,DiCabello2021,SupicBowles2020}.
Then the fact that any linear Bell inequality is maximized at some extremal point of the quantum set leads to a natural intuition that self-testing occurs on achieving maximal violation~\cite{SupicBowles2020,YN2013,Bamps-Pironio2015,Wang-Wu-Scarani2016,Coladangeloetal,Coopmanset.al.2019,Jed2020,Ishizaka2020,Rabelo2012,Baccarietal2020,DiCabello2021}.

In our work, we present a different approach for self-testing quantum states. By considering a Bell experiment called Hardy's test of nonlocality~\cite{Hardy1993} we show that self-test of a quantum state is possible also with a non-maximal violation. Such self-tests can be achieved for a two-parametric set of pure qubit states, which cover the whole spectrum from almost-not-entangled to almost-fully-entangled states. This result is interesting per se, as it allows self-testing of a broad spectrum of states through a single Bell test. On top of that, it has two important corollaries. First, the presented method using the concave cover approach is very general and can also be applied to other types of Bell tests. Second, as the behaviors which lead to self-testing are extremal points of the set of all quantum correlations, our result forms a stepping stone in the way of characterizing this set~\cite{Gohetal2018,Ishizaka2018,QVoid2019} and provides an important tool for future work in this area; one such step is realized, as we find that a part of the boundary of the quantum set is determined from Hardy's nonlocal correlations.

The subsequent parts of the paper are organized as follows. In Sec. II we first introduce Hardy’s test of nonlocality and then in Sec. III we characterize all two-qubit states that can demonstrate Hardy’s nonlocality. Sec. IV contain the derivation of the main results of this paper which follows by considering Hardy's test in a black-box scenario where the dimension of quantum state of the box is unknown. Finally, in the concluding Sec. V we provide a discussion and summary of our work.

\section{Hardy's Test}
The nonlocality test proposed by Hardy~\cite{Hardy1993} relates to a Bell experiment with two space-like separated parties, Alice and Bob, who share parts of a composite physical system. Alice randomly chooses to perform one of measurements $x\in\{A_0,A_1\}$ and Bob randomly chooses to perform one of measurements $y\in\{B_0,B_1\}$ on their respective parts.
Outcomes of all the measurements are binary, denoted $a\in\{\pm1\}$ for Alice and $b\in \{\pm 1\}$ for Bob.
Result of the experiment termed a behavior (correlation) is recorded in a vector of probabilities: $\vec{\mathscr{P}}=\{p(a,b\vert x,y): ~x\in \{A_0,A_1\},~y\in\{B_0,B_1\},~\mbox{and}~ a,b\in \{\pm1\}\}$. A behavior is termed \textit{local}, if it can be expressed in factorized form, {\it i.e.}, $p(a,b\vert x,y)=\int_{\lambda \in \Lambda} d\lambda~p(\lambda)p(a|x,\lambda)p(b|y,\lambda)$, where $p(\lambda)$ is a probability distribution over a set of local-hidden-variables $\Lambda$, and $p(a|x,\lambda),p(b|y,\lambda)$ are local response functions of Alice and Bob respectively, which without loss of any generality can be considered as deterministic \cite{Fine1982}.
Any behavior that cannot have a local-hidden-variable model is called \textit{nonlocal}.
Hardy showed that if the four conditions
\begin{subequations}
\begin{align}
	p_{\mbox{\tiny Hardy}}\equiv p(+1,~+1~\vert~ A_0,~B_0)&>0, \label{eq1a}\\
	p(+1,~-1~\vert~ A_0,~B_1)&=0, \label{eq1b}\\
	p(-1,~+1~\vert~ A_1,~B_0)&=0,\label{eq1c}\\
	p(+1,~+1~\vert~ A_1,~B_1)&=0,\label{eq1d}
\end{align}
\end{subequations}
are satisfied, then the resulting behavior is necessarily nonlocal. The probability $p_{\mbox{\tiny Hardy}}$ in Eq.~(\ref{eq1a}) quantifies the amount of nonlocality and it attains the maximum value {\small $p^{\tiny \mbox{max}}_{\tiny \mbox{Hardy}}=(-11 + 5 \sqrt{5})/2$} in quantum mechanics; the maximum is achieved with projective measurements on a pure two qubit state~\cite{Rabelo2012,SeshadreesanGhosh2011}.

\section{Two-qubit states showing Hardy's nonlocality}
In a two qubit state space {\small $\mathbb{C}^2\otimes \mathbb{C}^2$}, Hardy's nonlocal behaviors can result only from projective measurements on pure entangled states~\cite{GKar1997}. Due to local-unitary equivalence of measurements on suitably rotating the state by applying local unitary maps, without loss of generality,  let us consider projective measurements of the following form
\begin{subequations}
\begin{align}
A_0&=\vert0\rangle \langle 0\vert -\vert1\rangle \langle 1\vert,~~~~~A_1=\vert u_0\rangle\langle u_0\vert-\vert u_1\rangle\langle u_1\vert;\\
B_0&=\vert0\rangle \langle 0\vert -\vert1\rangle \langle 1\vert,~~~~~B_1~=\vert v_0\rangle\langle v_0\vert-\vert v_1\rangle \langle v_1\vert;
\end{align}\label{eq2}
\end{subequations}
where,
\begin{align*}
&\vert u_0\rangle=C_\alpha\vert 0\rangle + e^{\iota \phi}S_\alpha\vert 1\rangle,~\vert u_1\rangle=-S_\alpha\vert 0\rangle + e^{\iota \phi}C_\alpha\vert 1\rangle,\nonumber\\
&\vert v_0\rangle=C_\beta\vert 0\rangle + e^{\iota \xi}S_\beta\vert 1\rangle,~\vert v_1\rangle=-S_\beta\vert 0\rangle + e^{\iota \xi}C_\beta\vert 1\rangle, \nonumber
\end{align*}
with $C_z:=\cos(z/2),~S_z:=\sin(z/2)$ and $\alpha,\beta\in[0, \pi]$ and $\phi, \xi \in [0, 2\pi)$. Then, any pure two qubit state $\vert \psi\rangle$ satisfying constraints of Eqs.~(\ref{eq1b}-\ref{eq1d}) must be orthogonal to $\vert 0\rangle \!\otimes\!\vert v_1\!\rangle$, $\vert u_1\!\rangle \!\otimes\! \vert 0\rangle$, and $\vert u_0\!\rangle \!\otimes\! \vert v_0\!\rangle$, and therefore it must be of the form
\begin{align}
\vert \psi\rangle_{\tiny \mbox{Hardy}}\!=\!\frac{T_\alpha\vert u_0v_1\!\rangle\!+\! T_\beta\vert u_1v_0\!\rangle\!+\!\vert u_1v_1\!\rangle}{\sqrt{1+T^2_\alpha+T^2_\beta}};~~~T_z:=\tan\frac{z}{2}.\label{eq3}
\end{align}
\normalsize
The set of sixteen probabilities derived from the two qubit state in Eq. (\ref{eq3}) and measurements given by Eq.\eqref{eq2} can be expressed in an array as follows
\begin{widetext}
{ \begin{align}
\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}
\equiv
\begingroup
\setlength{\tabcolsep}{20pt}
\renewcommand{\arraystretch}{2}
\begin{array}{c||c|c|c|c|}
 & \!(\!+,+\!)\! & \!(\!+,-\!)\!&\!(\!-,+\!)\!&\!(\!-,-\!)\!\\ \hline\hline
\!A_0B_0\!\! &\!\dfrac{(\!1\!-\!r\!) r (\!1\!-\!s\!) s}{1\!-\!r s} \!&\!\dfrac{(\!1\!-\!r\!)^2 s}{1\!-\!r s} \!&\!(\!1\!-\!r\!)\!(\!1\!-\!s\!)\!&\!r \!\!\\\hline
\!A_0B_1\!\! & \!(1\!\!-\!\!r) s & 0 \!&\!\dfrac{(\!1\!-\!r\!) r s^2}{1\!-\!r s} \!&\!\dfrac{1\!-\!s}{1\!-\!r s}\!\!\\\hline
\!A_1B_0\!\! &\! \dfrac{(\!1\!-\!r\!) (\!1\!-\!s\!)}{1\!-\!r s} \!&\! \dfrac{r (\!1\!-\!s\!)^2}{1\!-\!r s} \!&\! 0 \!&\! s\!\!\\\hline
\!A_1B_1\!\! &\! 0 \!&\! 1-s \!&\!\dfrac{(\!1\!-\!r\!) s}{1\!-\!r s}\!&\! \dfrac{r (\!1\!-\!s\!) s}{1\!-\!r s}\!\!\\\hline
\end{array}
\endgroup
\label{eq4}
\end{align}}
\end{widetext}
where $r:=1-S^2_\alpha S^2_\beta$ and $s:=r^{-1}C^2_\alpha$.
Note that $0\leq r,s\leq 1$ for any choice of $\alpha,\beta \in [0,\pi]$, and a behavior $\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}$ is nonlocal if and only if $(r,s) \in (0,1)\times (0,1)$. The value of $\alpha$ and $\beta$ for a given value of $r$ and $s$ can be obtained from
{\small\begin{align}
\alpha=2\sin^{-1} \sqrt{1-rs},~~~~\beta=2\sin^{-1} \sqrt{(1-r)/(1-rs)}. \label{eq5}
\end{align}}
The two qubit state in Eq. (\ref{eq3}) when expressed in the standard basis, and in terms of the parameters $r$ and $s$, takes the following form
{\small	\begin{eqnarray}
\!\!\vert\!\psi(r,\!s)\!\rangle_{\tiny \mbox{Hardy}}\!\!=\!-\!\sqrt{\!\frac{(1\!-\!r) r (1\!-\!s) s}{1\!-\!r s}}~\vert00\rangle\!&-& e^{\iota \xi}\!\sqrt{\!\frac{(1\!-\!r)^2 s}{1\!-\!r s}}~\vert 01\rangle \nonumber \\
	-e^{\iota \phi}\sqrt{\!(\!1\!-\!r\!) (\!1\!-\!s\!)}~\vert 10\rangle\!&+&e^{\iota (\xi+\phi)}\sqrt{\!r}~\vert 11\rangle.\label{eq6}
\end{eqnarray}}
Hardy nonlocal states in Eq.(\ref{eq3}) [or Eq.(\ref{eq6})] covers all (up to local-unitary rotations of basis) pure two qubit entangled states except the maximally entangled ones~\cite{Goldstein1994,Jordan1994}. Moreover, in the two qubit space, due to the constraints in Eq.~(\ref{eq1a}-\ref{eq1b}), these states are uniquely determined from any arbitrarily fixed measurements~\cite{Jordan1994}.

\section{Hardy’s test with unknown state of arbitrary dimension}
Let us now consider a black-box experiment (under i.i.d. assumption) where the quantum state and measurements of the box are unknown to Alice and Bob. The main result of this paper is that if the probabilites from the black-box are in the form given by Eq.~(\ref{eq4}), then they are a self-test of the two qubit state in Eq.~(\ref{eq6}). The black-box experiment records $p(-1,-1\vert A_0,B_0)\equiv\mathbf{r}$ and $p(-1,-1\vert A_1,B_0)\equiv\mathbf{s}$, and then verifies if the remaining probabilities are expressible in terms of $\mathbf{r}$ and $\mathbf{s}$ in the form of Eq.(\ref{eq4}); if true, then the claim is that, it is a self-test of the two qubit state $\vert \psi(\mathbf{r},\mathbf{s})\rangle_{\tiny \mbox{Hardy}}$ in Eq.~(\ref{eq6}).

The main idea in the proof of our claim follows from a Jordan canonical form for measurements $A_0, A_1$ ($B_0,B_1$) of Alice (Bob), and application of Jensen's inequality. Suppose some unknown state $\rho$ is shared between Alice and Bob and $\Pi_{a\vert x}$ ($\Pi_{b\vert y}$) is the measurement operator associated with outcome $a$ ($b$) when Alice (Bob) measures observable $x$ ($y$).
Then we have, $p(a,b\vert x,y)=\mbox{Tr}(\rho~ \Pi_{a\vert x} \otimes \Pi_{b\vert y})$.
Since the dimension of the state space is unrestricted, Neumark's dilation theorem allows us to reduce the analysis to projective measurements.
Then, let the observable of Alice and Bob be some Hermitian operators with eigenvalues in $\{+1,-1\}$ as follows
{\small \begin{eqnarray*}
	x&=& (+1)~\Pi_{+\vert x} ~+~(-1)~\Pi_{-\vert x}~~~~\mbox{where}~x\in\{A_0,A_1\},\\
	y&=& (+1)~\Pi_{+\vert y} ~+~(-1)~\Pi_{-\vert y}~~~~\mbox{where}~y\in\{B_0,B_1\}.
\end{eqnarray*}}
These observables can be written in a Jordan canonical form derived in Ref.~\cite{Masanes2006} (also see~\cite{Masanes2005,Pironio2009}), which states that: for any set of four projection operators {\small $ \{\Pi_{+\vert \mathcal{M}_0},~\Pi_{-\vert \mathcal{M}_0},~\Pi_{+\vert \mathcal{M}_1},~\Pi_{-\vert \mathcal{M}_1}\}$} acting on a Hilbert space $\mathcal{H}$ and satisfying conditions {\small $\Pi_{+\vert \mathcal{M}_0} +\Pi_{-\vert \mathcal{M}_0}=I$} and {\small $\Pi_{+\vert \mathcal{M}_1} +\Pi_{-\vert \mathcal{M}_1}=I$}, there is an orthonormal basis of {\small $\mathcal{H}$} in which all the four operators are simultaneously block diagonal with maximum block size {\small $2\times 2$}. The orthonormal basis induces a direct sum decomposition {\small $\mathcal{H}=\oplus_k\mathcal{H}^k$} where dimension of each component subspace {\small $\mathcal{H}^k$} is at most two.
Then, the four projection operators can be decomposed as {\small$\Pi_{\pm\vert \mathcal{M}_{0(1)}}=\oplus_k~\Pi_{\pm\vert \mathcal{M}_{0(1)}}^k$} and each component {\small$\Pi_{\pm\vert \mathcal{M}_{0(1)}}^k$} acts on the subspace $\mathcal{H}^k$.
The projector on subspace $\mathcal{H}^k$ can be written as {\small$\Pi^k=\Pi_{+\vert \mathcal{M}_{0}}^k+\Pi_{-\vert \mathcal{M}_{0}}^k=\Pi_{+\vert \mathcal{M}_{1}}^k+\Pi_{-\vert \mathcal{M}_{1}}^k$}.
On applying the stated result to the observables {\small$A_0,A_1~(B_0,B_1)$} and state space {\small$\mathcal{H}_A ~\left(\mathcal{H}_B\right)$} of Alice (Bob) gives
{\small \begin{equation}
p(\!a,\!b\vert x,\!y\!)=\!\sum_{i,j} \mu_{ij}\mbox{Tr}(\rho_{ij}\Pi^i_{a\vert x}\!\! \otimes\! \Pi^j_{b\vert y})\equiv \!\sum_{i,j} \mu_{ij}p_{ij}(\!a,\!b\vert x,\!y\!), \label{eq7}
\end{equation}}
where {\small $\mu_{ij}\!=\!\mbox{Tr}(\rho\Pi^i\!\otimes\!\Pi^j)$} and satisfies {\small$\sum_{i,j}\mu_{ij}=1$} and {\small $\mu_{ij}\geq 0$}, and {\small$\rho_{ij}\!=\!(\Pi^i\!\otimes\!\Pi^j\rho\Pi^i\!\otimes\!\Pi^j)/\mu_{ij}$} is trace one positive operator on component subspace {\small $\mathcal{H}_A^i\otimes \mathcal{H}_B^j$}.

Now let us define from the probabilities in Eq.\eqref{eq4} a function {\small $\Omega(\mathbf{r},\mathbf{s})=\sum_{a,b,x,y} c_{abxy} ~p(a,b\vert x,y)+c_0$} where $c_{abxy}$ and $c_0$ are some real coefficients. On applying Eq.(\ref{eq7}) to each probability term in {\small $\Omega(\mathbf{r},\mathbf{s})$} we get
{\small\begin{align}
\Omega (\mathbf{r},\mathbf{s})= \sum_{i,j} \mu_{ij}~ \Omega (r_{ij}, ~s_{ij}), \label{eq8}
\end{align}}
where {\small $r_{ij}\!=\!p_{ij}(\!-1,\!-1\vert A_0,B_0\!)$} and {\small $s_{ij}\!=\!p_{ij}(\!-1,\!-1\vert A_1,B_0\!)$}. From Eq.(\ref{eq7}) one can find that $\mathbf{r}=\sum_{i,j} \mu_{ij}r_{ij}$ and $\mathbf{s}=\sum_{i,j} \mu_{ij}s_{ij}$. Furthermore, when the black-box statistics satisfies the zero constraints of Hardy's test, then the same is true in every {\small $\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace. Now we like to state and prove the following Lemma.
\begin{lem}
 Let {\small$\mathcal{E}(\mathbf{r},\!\mathbf{s})\!: \!(0,\!1)\!\times\!(0,\!1)\!\rightarrow \!\mathbb{R}$} be a concave cover of {\small$\Omega (\mathbf{r},\mathbf{s})$}, and suppose $\mathcal{R}$ is the set of points from the domain for which {\small$\mathcal{E}(\mathbf{r},\mathbf{s})=\Omega (\mathbf{r},\mathbf{s})$}. Then for all {\small$(\mathbf{r},\mathbf{s})\in \mathcal{R}$}, and for all {\small$(i,j)$},  {\small$r_{ij}=\mathbf{r}$} and {\small$s_{ij}=\mathbf{s}$}, provided {\small$\Omega(\mathbf{r},\mathbf{s})$} is a strictly concave function of {\small$\mathbf{r}$} and {\small$\mathbf{s}$} over the region {\small$\mathcal{R}$}. \label{Lemma1}
\end{lem}
\begin{proof}
A \textit{concave cover} {\small$\mathcal{E}(\mathbf{r},\mathbf{s})$} for function {\small$\Omega(\mathbf{r},\mathbf{s})$} over its domain {\small$(0,1)\times(0,1)$} is defined as the lowest-valued concave function that overestimates or equals $\Omega(\mathbf{r},\mathbf{s})$ in its domain.
Now consider Jensen's inequality which states that for any concave real function {\small$f(x):\mathbb{R}^n\rightarrow \mathbb{R},~f(\sum_{k=1}^m p_k x_k) \geq \sum_{k=1}^m p_k f(x_k)$} where {\small$x_k\in\mathbb{R}^n$} and {\small$p_k\geq 0$} for all {\small$k$}, and {\small$\sum_kp_k=1$}. When {\small$f(x)$} is nonlinear, equality holds if and only if {\small$x_1=x_2=...=x_m$}. On applying Jensen's inequality to the function {\small $\mathcal{E}(\mathbf{r},\mathbf{s}):\mathbb{R}^2\rightarrow \mathbb{R}$} we get
{\small\begin{eqnarray}
\mathcal{E}(\mathbf{r},\mathbf{s})=\mathcal{E}(\sum_{i,j} \mu_{ij}~( r_{ij}, ~s_{ij}))
\geq \sum_{ij} \mu_{ij}~ \mathcal{E}(r_{ij},~ s_{ij}). \label{eq9}
\end{eqnarray}}
On the other hand, when {\small $(\mathbf{r},\mathbf{s}) \in \mathcal{R}$} we have {\small$\mathcal{E}(\mathbf{r},\mathbf{s}) =\Omega (\mathbf{r},\mathbf{s})$}, and since {\small$ \Omega(r_{ij},s_{ij})\leq\mathcal{E}(r_{ij},s_{ij})$}, using Eq.\eqref{eq10}, we get
{\small\begin{align}
\mathcal{E}(\mathbf{r},\mathbf{s}) =\sum_{i,j} \mu_{ij}~\Omega(r_{ij}, s_{ij})
 \leq \sum_{i,j} \mu_{ij}~\mathcal{E}(r_{ij}, s_{ij}).\label{eq10}
\end{align}}
Then, Eq.~(\ref{eq9}) and Eq.~(\ref{eq10}) imply that if $(\mathbf{r},\mathbf{s})\in \mathcal{R}$
{\small\begin{align}
\mathcal{E}(\sum_{i,j} \mu_{ij}( r_{ij}, ~s_{ij})) = \sum_{ij} \mu_{ij}~ \mathcal{E}(r_{ij}, s_{ij}). \label{eq11}
\end{align}}
Thus in Jensen's inequality~(\ref{eq9}), equality holds in the region $\mathcal{R}$. Therefore, if $\Omega(\mathbf{r},\mathbf{s})$ is a strictly concave function of $\mathbf{r}$ and $\mathbf{s}$ over the region $\mathcal{R}$, in every {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace value of the parameters  $(r_{ij},s_{ij})$ must be same, i.e., $r_{ij}=\mathbf{r}$ and $s_{ij}=\mathbf{s}$ for all $i,j$.
\end{proof}
\begin{figure}[t!]
	\begin{center}
		\includegraphics[angle=0, width=0.3\textwidth]{Fig1.jpg}
		\caption{ In the light (blue) shaded region, {\small$\Omega^{\ast}(\mathbf{r},\mathbf{s})$} is concave. {\small$\mathcal{R}^{\ast}(\mathbf{r},\mathbf{s})$} is the dark (blue) shaded region,  where the concave cover {\small$\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})=\Omega^{\ast}(\mathbf{r},\mathbf{s})$. } The dark (blue) dot shows point $(\mathbf{r},\mathbf{s})=(\frac{\sqrt{5}-1}{2} ,\frac{\sqrt{5}-1}{2})$ where $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ is maximum.}  \label{Fig1}
	\end{center}
\end{figure}
Let us now give an example where Lemma \ref{Lemma1} is applicable. Consider a function defined by the probability which quantifies success in Hardy's test, i.e.,
{\small\begin{align}
\Omega^{\ast}(\mathbf{r},\mathbf{s})\equiv p(+1,+1\vert A_0,B_0)=\frac{\mathbf{r}(1\!-\!\mathbf{r})\mathbf{s}(1\!-\!\mathbf{s})}{1\!-\!\mathbf{r} \mathbf{s}}.
\label{eq12}
\end{align}}
The function {\small$\Omega^{\ast}(\mathbf{r},\mathbf{s})$} is a nonlinear function and it is concave in a part of its domain. Further, for the considered function, there exist a concave cover {\small $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})$}, and a region {\small $\mathcal{R}^{\ast}\subset (0,1)\times(0,1)$} where {\small$\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})=\Omega^{\ast}(\mathbf{r},\mathbf{s})$}.
We constructed the cover function numerically and then find the region $\mathcal{R}^{\ast}$ as shown in Fig. \ref{Fig1} (see Appendix~(A) for the details of the computational method).
From the considered example and application of Lemma~\ref{Lemma1}, we conclude that if probability distribution of the black-box is such that $(\mathbf{r},\mathbf{s})\in \mathcal{R}^{\ast}$ then in every subspace  {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$, $r_{ij}=\mathbf{r}$ and $s_{ij}=\mathbf{s}$}. We now ask whether such a property can hold for all $(\mathbf{r},\mathbf{s})\in (0,1)\times (0,1)$. In other words, is it possible to vary over the possible choice of functions $\Omega(\mathbf{r},\mathbf{s})$ such that union of all the resulting regions $\mathcal{R}$ covers the full parameter space {\small$(0,1)\times (0,1)$}? Such a property will free the black box parameters $\mathbf{r}$ and $\mathbf{s}$ from any restrictions and lead to an interesting extension of Lemma~\ref{Lemma1}. We find answer to the question in affirmative through the following lemma.
\begin{lem}
 If the black-box statistics in the Bell-experiment is of the form $\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}(\mathbf{r},\mathbf{s})$ as in Eq.~(\ref{eq4}), then for all $(\mathbf{r},\mathbf{s})\in(0,1)\times (0,1)$ and in all $\mathcal{H}_A^i\otimes \mathcal{H}_B^j$ subspace, $r_{ij}=\mathbf{r}$ and $s_{ij}=\mathbf{s}$. \label{Lemma2}
 \end{lem}
\begin{proof}
First we recall that Lemma \ref{Lemma1} is applicable to any function {\small $\Omega(\mathbf{r},\mathbf{s})$} defined as some linear combination of all the probabilities in Eq~(\ref{eq4}) plus a constant term. We find that, a proof of Lemma~\ref{Lemma2} follows on considering simply a single parameter family of functions
{\small \begin{align}
\Omega_{\nu}(\mathbf{r},\!\mathbf{s}) \!&=\!\Omega^{\ast}(\mathbf{r},\!\mathbf{s}) \!+\nu~\!p(\!+\vert A_0\!) \!+(1\!-\!\nu)~\!p(\!-\vert A_0\!)\!-\!1/2,\nonumber\\
	\!&=\!\Omega^{\ast}(\mathbf{r},\!\mathbf{s})\!+\nu~\!(\mathbf{s}\!-\!\mathbf{r}\mathbf{s}) \!+(1\!-\!\nu)~\!(1\!-\!\mathbf{s}\!+\!\mathbf{r}\mathbf{s})\!-\!1/2, \label{eq13}
\end{align}}
where {\small $0\leq \nu \leq 1$}, {\small $p(\!+\vert \!A_0\!)\!=\!p(\!+,\!+\vert A_0,\!B_0\!)+p(\!+,\!-\vert A_0,\!B_0\!)$, and $p(\!-\vert A_0\!)=1-p(\!+\vert A_0\!)$}. Note that at {\small$\nu =1/2$}, {\small$\Omega_{1/2}(\mathbf{r},\mathbf{s})\equiv \Omega^{\ast}(\mathbf{r},\mathbf{s})$} and {\small $\mathcal{R}_{1/2}\equiv \mathcal{R}^{\ast}$}.
For the family of functions {\small$\Omega_{\nu}(\mathbf{r},\mathbf{s})$} on considering {\small$\nu \in \{d/N:~d\in\{1,2,...,N-1\} \}$}, we find that for a sufficiently large value of {\small$N$}, {\small$\cup_{\nu} \mathcal{R}_{\nu}= (0,1)\times (0,1)$}. To observe the proof, first let us see  Fig.~\ref{Fig2} where the region {\small$\cup_{\nu} \mathcal{R}_{\nu}$} is shown when for {\small$N=10$}, i.e., when the parameter $\nu$ takes values from the set $\{0.1,~0.2,...,0.9\}$. Note that the region where {\small$r_{ij}=\mathbf{r}$} and {\small$s_{ij}=\mathbf{s}$}, in all {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace, is expanded. A more detailed exposition of the proof which follows on increasing the value of {\small$N$} is provided in the Appendix~(B).
\end{proof}
Finally, in the following we prove the main result of the paper by using Lemma~\ref{Lemma2}.
\begin{figure}[t!]
	\begin{center}
		\includegraphics[angle=0, width=0.3\textwidth]{Fig2.jpg}
		\caption{ The dark (blue) shaded region {\small$\cup_{\nu} \mathcal{R}_{\nu}$} is shown when the parameter {\small$\nu$} takes values from the set {\small$\{0.1,~0.2,...,~0.9\}$}. If {\small$(\mathbf{r},\mathbf{s})$} belongs to the shaded region, then {\small$(r_{ij},s_{ij})=(\mathbf{r},\mathbf{s})$} in every {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace. When {\small$\cup_{\nu} \mathcal{R}_{\nu}$} is derived from {\small$\{\nu=k/N:~k\in\{1,2,...,N-1\}\}$}, as {\small$N$} becomes sufficiently large the shaded region covers all {\small$(\mathbf{r},\mathbf{s}) \in (0,1)\times (0,1)$}.}  \label{Fig2}
	\end{center}
\end{figure}
\begin{thm*}
In a black-box Bell experiment, if a behavior $\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}(\mathbf{r},\mathbf{s})$ of the form given in Eq.~(\ref{eq4}) is observed, then the state of unknown system $\rho^{AB}$ is equivalent up to local isometries to $ \zeta^{AB}\otimes\vert\psi(\mathbf{r},\mathbf{s})\rangle_{\tiny \mbox{Hardy}}^{A'B'}\langle\psi(\mathbf{r},\mathbf{s})\vert$, where $\vert\psi(\mathbf{r},\mathbf{s})\rangle^{A'B'}_{\tiny \mbox{Hardy}}$ is the pure two qubit Hardy state given by the Eq.~(\ref{eq6}) and $\zeta^{AB}$ is some arbitrary bipartite state.
\end{thm*}
\begin{proof}
Consider a purification $\vert\chi\rangle^{ABP}$ of the unknown state $\rho^{AB}$, in brief we denote such a purification simply as $\vert\chi\rangle^{AB}$ since the desired local isometry is not required to act on the purification space~\cite{SupicBowles2020}. Also observables $A_0,A_1,B_0,B_1$ can be considered in a $2\times2$ block diagonal form~\cite{SupicBowles2020}. Then, there is a basis in which the observables are in the following block diagonal form
{\small\begin{align}
\Pi^i_{+\vert A_0}&=\vert2i\rangle\langle2i\vert,~~~~~~\Pi^i_{-\vert A_0}=\vert2i+1\rangle\langle2i+1\vert;\nonumber\\
\Pi^i_{+\vert A_1}&=\vert u_{2i}\rangle\langle u_{2i}\vert,~~~
\Pi^i_{-\vert A_1}=\vert u_{2i+1}\rangle\langle u_{2i+1}\vert;\nonumber\\
\Pi^j_{+\vert B_0}&=\vert2j\rangle\langle2j\vert,~~~~~~ \Pi^j_{-\vert B_0}=\vert2j+1\rangle\langle2j+1\vert;\nonumber\\
\Pi^j_{+\vert B_1}&=\vert v_{2j}\rangle\langle v_{2j}\vert,~~~ \Pi^j_{-\vert B_1}=\vert v_{2j+1}\rangle\langle v_{2j+1}\vert;
\label{eq14}
\end{align}}
where,
{\small \begin{align}
\vert u_{2i}\!\rangle\!&=\!C_{\alpha_i}\vert 2i\rangle \!+ \!e^{\mathtt{i} \phi}S_{\alpha_i}\vert 2i\!+\!1\rangle,~
\vert u_{2i+1}\!\rangle\!=\!-S_{\alpha_i}\vert 2i\rangle \!+\! e^{\mathtt{i} \phi}C_{\alpha_i}\vert 2i\!+\!1\rangle; \nonumber\\
\vert v_{2j}\!\rangle\!&=\!C_{\beta_j}\vert 2j\rangle \!+\! e^{\mathtt{i} \xi}S_{\beta_j}\vert 2j\!+\!1\rangle,~
\vert v_{2j+1}\!\rangle\!=\!-S_{\beta_j}\vert 2j\rangle \!+\! e^{\mathtt{i} \xi}C_{\beta_j}\vert 2j\!+\!1\rangle, \nonumber
\end{align}}
with $i,j\in \{0,1,2,...\}$, $\alpha_i,\beta_j\in[0, \pi]$ and $\phi, \xi \in [0, 2\pi)$.
Note that, without loss of generality, the phases $\xi$ and $\phi$ can be considered independent of indices $i$ and $j$ in every {\small $\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace since this can be achieved by choosing a suitable measurement basis and rotating the state by local unitaries such that the probabilities $p_{ij}(a,b\vert x,y)$ remain invariant. Also note that in the subspace {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$ parameters $r_{ij}$ and $s_{ij}$} are given by {\small$r_{ij}=1-S^2_{\alpha_i}S^2_{\beta_j},~~s_{ij}=r^{-1}_{ij}C^2_{\alpha_i}$}.
Now suppose $\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}(\mathbf{r},\mathbf{s})$ is observed in the black-box experiment, then Lemma~\ref{Lemma2} implies that in each subspace {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j $}, $r_{ij}=\mathbf{r}$ and $s_{ij}=\mathbf{s}$ and, therefore, $\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}(r_{ij},s_{ij})=\vec{\mathscr{P}}_{\tiny \mbox{Hardy}}(\mathbf{r},\mathbf{s})$. This can be true if and only if in {\small$\mathcal{H}_A^i\otimes \mathcal{H}_B^j$} subspace, $\rho_{ij}=\vert\psi_{\tiny \mbox{Hardy}}\rangle_{ij}\langle \psi_{\tiny \mbox{Hardy}}\vert$, where
{\small \begin{align}
&\vert\!\psi_{\tiny \mbox{Hardy}}\!\rangle_{ij}\!=\!-\sqrt{\!\frac{(1\!-\!\mathbf{r}) \mathbf{r} (1\!-\!\mathbf{s}) \mathbf{s}}{1\!-\!\mathbf{r} \mathbf{s}}}~\vert2i,2j\rangle\!-\!e^{\iota\xi}\!\sqrt{\!\frac{(1\!-\!\mathbf{r})^2 \mathbf{s}}{1\!-\!\mathbf{r} \mathbf{s}}}~\vert 2i,2j\!+\!1\rangle \nonumber\\
&-e^{\iota\phi}\!\sqrt{\!(1\!-\!\mathbf{r}) (1\!-\!\mathbf{s})}~\vert 2i\!+\!1,2j\rangle\!+\!e^{\iota(\xi+\phi)}\!\sqrt{\!\mathbf{r}}~\vert 2i\!+\!1,2j\!+\!1\rangle,
\label{eq15}
\end{align}}
Hence, the unknown state {\small$\vert\chi \rangle^{AB}$} can only be a direct sum of the form {\small$\vert\chi \rangle^{AB}=\bigoplus_{i,j}\sqrt{\mu_{ij}}~\vert\psi_{\tiny \mbox{Hardy}}\rangle_{ij}$}. Finally, we can give local isometries $\Phi^A$ and $\Phi^B$ such that
{\small\begin{equation*}
	(\Phi^A\otimes\Phi^B)\vert\chi\rangle^{AB}\vert 00\rangle^{A'B'}=\zeta^{AB}\otimes\vert\psi(\mathbf{r},\mathbf{s})\rangle_{\tiny \mbox{Hardy}}^{A'B'}\langle\psi(\mathbf{r},\mathbf{s})\vert,
\end{equation*}}
where components of the $\vert 00\rangle^{A'B'}$ are local ancilla qubits appended to the unknown state $\vert\chi\rangle^{AB}$, and after application of the local isometry $\Phi^A\otimes\Phi^B$ we want to get the target state $\vert\psi(\mathbf{r},\mathbf{s})\rangle_{\tiny \mbox{Hardy}}^{A'B'}$ along with some arbitrary bipartite state $ \zeta^{AB}$.
The isometry map $\Phi^A=\Phi^B=\Phi$ with $\Phi~\vert 2k,0\rangle^{XX'} \mapsto \vert 2k,0\rangle^{XX'}$ and $ \Phi~\vert 2k+1,0\rangle^{XX'} \mapsto \vert 2k,1\rangle^{XX'}$, where $XX'\in\{AA',BB'\}$, has the desired property.
This concludes our proof.
\end{proof}
The Theorem implies that any behavior in the form of Eq.~(\ref{eq4}) are extremal points of the set of quantum behaviors and hence it determines a part of the quantum boundary. Further, one can check that all the behaviors that are Hardy nonlocal but not in the form given by Eq.~(\ref{eq4}) are in the interior of the quantum set and they must arise on measuring mixed entangled states in higher than two qubit dimensions.

\section{Conclusion}
Hardy's test of nonlocality, viewed differently, is a form of tailored Bell inequality~\cite{Acin2012, Alexia2017}. However, in distinction to, for instance, \textit{tilted}-Clauser-Horne-Shimony-Holt inequality~\cite{Acin2012}, Hardy's test places certain constraints on some outcome probabilities. Importance of tailored Bell inequalities is starkly revealed through various protocols for self-testing, randomness generation, quantum key distribution, etc., wherever the knowledge of the geometry of the quantum set of correlations plays a crucial role~\cite{Gohetal2018}. In this sense tailoring by constraining outcome probabilities can possibly play an important role, as we have demonstrated in this work through the Hardy's test. Here we note that for self-testing two qubit states and measurements on them, tilted-CHSH inequality is more powerful than Hardy's test, however, Hardy's correlations can self-test almost all two qubit states from the whole spectrum of violations, ranging from arbitrarily small to maximal violation in the Hardy's test of nonlocality. The fact that our approach can have broader applications can be seen by noticing the three key features in the derivation of our results: (i) due to the constraints on probabilities a Hardy nonlocal state, with local state space $\mathbb{C}^2$, is pure and unique for an arbitrarily fixed measurement~\cite{GKar1997,Jordan1994}, (ii) existence of a simple Jordan canonical form for two dichotomic observables for each party~\cite{Masanes2006,Masanes2005,Pironio2009}, and (iii) possibility of finding a concave cover to an arbitrary linear combination of outcome probabilities and application of Jensen's inequality. Thus, there can be tailored Bell tests other than Hardy's where the three properties may hold. In general, these three features may be found in any Bell scenario with $n$-parties, two measurements per party, and two outcomes to every measurement.

To summarize, in this work, we have presented a method for self-testing quantum states by using the concave cover of a linear combination of observed outcome probabilities. We showcased the applicability of this approach on Hardy's test of nonlocality, leading to a two-parametric set of characterizable states with a broad spectrum of entanglement properties. This technique allows a full (up to local isometries) specification of the source state, even without maximal violation in Hardy's test. Our results show many potential development paths as the method introduced in this work can be possibly applied to different kind of Bell tests, allowing better specification of source states in a device-independent fashion. At the same time, it provides a very powerful tool for characterizing the boundary of the quantum set of correlations.

\begin{acknowledgments}
AR is supported by the National Research Foundation of Korea (NRF-2021R1A2C2006309), an Institute of Information and Communications Technology Promotion (IITP) grant funded by the Korean government (MSIP) (Grant No. 2019-0-00831) and the University IT Research Center (ITRC) Program (IITP-2021-2018-0-01402). AR, M. Pivoluska, and M. Plesch acknowledge funding and support from VEGA Project~No.~2/0136/19. M. Pivoluska, and M. Plesch additionally acknowledge GAMU project MUNI/G/1596/2019.  M.B. acknowledges funding from the
National Mission in Interdisciplinary Cyber-Physical systems from the Department of Science and Technology through
the I-HUB Quantum Technology Foundation (Grant No. I-HUB/PDF/2021-22/008), support through the research grant of INSPIRE Faculty fellowship from the Department of Science and Technology, Government of India, and the start-up research grant from SERB, Department of Science and Technology (Grant No. SRG/2021/000267).

AR thankfully acknowledge fruitful discussions and feedback from  Prof. M.J.W. Hall, Prof. Joonwoo Bae, Prof. Guruprasad Kar, and Dr. Ramij Rahaman at different stages of this work.

\end{acknowledgments}

\begin{widetext}

\section*{Appendix~(A): Computation of Concave cover for functions $\boldmath{\Omega(\mathbf{r},\mathbf{s})}$} \label{apnA}
Let us describe our method for computing the concave covers to the class of functions $\Omega(\mathbf{r},\mathbf{s})$ introduced in  Lemma-1 of the main text. Consider the first concrete example to which application of the Lemma-1 gives interesting insight,
\begin{align}
	\Omega^{\ast}(\mathbf{r},\mathbf{s})=\frac{\mathbf{r}(1\!-\!\mathbf{r})\mathbf{s}(1\!-\!\mathbf{s})}{1\!-\!\mathbf{r} \mathbf{s}}.
\end{align}
The three dimensional plot of $	\Omega^{\ast}(\mathbf{r},\mathbf{s})$ is shown in Fig.(\ref{SupFig1a}) which indicates that the function can be concave in some parts of its domain. On computing the two eigenvalues of Hessian matrix of $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ we find that there a region of domain where both the eigenvalues are negative, therefore, the function is concave in this region; this region is shown in Fig.(\ref{SupFig1b}).

\begin{figure}[!ht]
	\begin{tabularx}{\linewidth}{CC}
		\begin{subfigure}[b]{0.4\textwidth}
			\includegraphics[width=\linewidth]{Fig3a.jpg}
			\caption{A $3$-dimensional  plot of $\Omega^{\ast}(\mathbf{r},\mathbf{s})$}
			\label{SupFig1a}
		\end{subfigure}
		&
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\linewidth]{Fig3b.jpg}
			\caption{In shaded (light blue) region $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ is concave; bold (blue) dot in the middle is the point $(\mathbf{r},\mathbf{s})=(\frac{\sqrt{5}-1}{2} ,\frac{\sqrt{5}-1}{2})$ where the function takes the maximum value.}
			\label{SupFig1b}
		\end{subfigure}

	\end{tabularx}
	\caption{Region of concavity for the function  $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ }
	\label{SupFig1}
\end{figure}
Next, we constructed a concave cover $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})$ of the function $\Omega^{\ast}(\mathbf{r},\mathbf{s})$. Recall that by definition a concave cover of the function $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ is the lowest-valued concave function that overestimates or equals $\Omega^{\ast}(\mathbf{r},\mathbf{s})$ in its domain. We numerically computed the concave cover by applying the ConvexHullMesh function provided in MATHEMATICA. The concave cover $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})$ is shown in Fig.(\ref{SupFig2a}) and the part of domain $\mathcal{R}^{\ast}$ where $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})=\Omega^{\ast}(\mathbf{r},\mathbf{s})$ is shown as a shaded region in Fig.(\ref{SupFig2b}). Figure one of the main text is obtained by superposing Fig.(\ref{SupFig1b}) and Fig.(\ref{SupFig2b}).

By a similar procedure one can construct the concave cover $\mathcal{E}(\mathbf{r},\mathbf{s})$ to any function $\Omega(\mathbf{r},\mathbf{s})$, and then the corresponding region $\mathcal{R}$ where $\mathcal{E}(\mathbf{r},\mathbf{s})=\Omega(\mathbf{r},\mathbf{s})$.
\begin{figure}[t!]
	\begin{tabularx}{\linewidth}{CC}
		\begin{subfigure}[b]{0.4\textwidth}
			\includegraphics[width=\linewidth]{Fig4a.jpg}
			\caption{Concave cover $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})$ of function $\Omega^{\ast}(\mathbf{r},\mathbf{s})$}
			\label{SupFig2a}
		\end{subfigure}
		&
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\linewidth]{Fig4b.jpg}
			\caption{Shaded (Blue) region shows $\mathcal{R}^{\ast}$ where $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})=\Omega^{\ast}(\mathbf{r},\mathbf{s})$}
			\label{SupFig2b}
		\end{subfigure}

	\end{tabularx}
	\caption{Construction of concave cover $\mathcal{E}^{\ast}(\mathbf{r},\mathbf{s})$ and region $\mathcal{R}^{\ast}$}
	\label{SupFig2}
\end{figure}

\section*{ Appendix~(B): Details of the proof for Lemma-$2$ in the main text}  \label{apnB}

Let us now consider the following single parameter family of functions, of the two variables $\mathbf{r}$ and $\mathbf{s}$, introduced in the construction of the proof of Lemma-2 in the main text
\begin{align}
	\Omega_{\nu}(\mathbf{r},\!\mathbf{s}) \!&=\!\Omega^{\ast}(\mathbf{r},\!\mathbf{s}) \!+\nu~\!p(\!+1\vert A_0\!) \!+(1\!-\!\nu)~\!p(\!-1\vert A_0\!)\!-\!1/2,\nonumber\\
	\!&=\!\Omega^{\ast}(\mathbf{r},\!\mathbf{s})\!+\nu~\!(\mathbf{s}\!-\!\mathbf{r}\mathbf{s}) \!+(1\!-\!\nu)~\!(1\!-\!\mathbf{s}\!+\!\mathbf{r}\mathbf{s})\!-\!1/2, \label{Eq1}
\end{align}
where the parameter $\nu\in (0,1)$ and variables $(\mathbf{r},~\mathbf{s})\in (0,1)\times(0,1)$. Firstly, we constructed finite subset of functions from the family in~Eq.(\ref{Eq1}) containing $N \!- \!1$ members: $S_{N}=\{\Omega_{\nu}(\mathbf{r},\mathbf{s}):~\nu \in \{d/N:~d\in\{1,2,...,N \!-\! 1 \} \}\}$. Let us denote the $k$-th member of $S_{N}$ by $S_{N}(k)$ where $k\in\{1,2,,...,N \!-\! 1\}$.
\begin{figure}[h!]
	\begin{tabularx}{\linewidth}{CCCC}
		\centering
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$\mathcal{R}_{\frac{1}{4}}$~$(\nu=1/4)$}
			\includegraphics[width=\textwidth]{Fig5a.jpg}
			\label{Fig1a}
		\end{subfigure}

		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$\mathcal{R}_{\frac{1}{2}}$~$(\nu=1/2)$}
			\includegraphics[width=\textwidth]{Fig5b.jpg}
			\label{Fig1b}
		\end{subfigure}

		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$\mathcal{R}_{\frac{3}{4}}$~$(\nu=3/4)$}
			\includegraphics[width=\textwidth]{Fig5c.jpg}
			\label{Fig1c}
		\end{subfigure}

		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$\cup_{\nu} \mathcal{R}_{\nu}$}
			\includegraphics[width=\textwidth]{Fig5d.jpg}
			\label{Fig1d}
		\end{subfigure}

	\end{tabularx}
	\caption{Plots of regions $\mathcal{R}_{\nu}$ and $\cup_{\nu}\mathcal{R}_{\nu}$ for $\nu\in\{\frac{1}{4},\frac{1}{2},\frac{3}{4}\}$ and $N=4$.}
	\label{SupFig3}
\end{figure}

\begin{figure}[h!]
	\begin{tabularx}{\linewidth}{CCCC}
		\centering
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$N=2$}
			\includegraphics[width=\textwidth]{Fig6a.jpg}
			\label{Fig2a}
		\end{subfigure}
		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$N=5$}
			\includegraphics[width=\textwidth]{Fig6b.jpg}
			\label{Fig2b}
		\end{subfigure}

		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$N=10$}
			\includegraphics[width=\textwidth]{Fig6c.jpg}
			\label{Fig2c}
		\end{subfigure}

		&
		\begin{subfigure}[b]{0.25\textwidth}
			\centering
			\caption{$N=100$}
			\includegraphics[width=\textwidth]{Fig6d.jpg}
			\label{Fig2d}
		\end{subfigure}

	\end{tabularx}
	\caption{Plots of region $\cup_{\nu}\mathcal{R}_{\nu}$ for different values of $N$.}
	\label{SupFig4}
\end{figure}
Secondly, we numerically constructed concave envelope functions of $S_{N}(k)$ for all values of $k$ by ConvexHullMesh function
in MATHEMATICA, then the region of domain where value of a concave envelope function coincide with the value of $S_{N}(k)$ was found as $\mathcal{R}_{\nu}$. Finally, for fixed values of $N$ the region $\cup_{\nu} \mathcal{R}_{\nu}$ was constructed. For example, the Fig.(\ref{SupFig3}) shows these region for $N=4$. By the above procedure, for a sufficiently large value of $N$ we find that $\cup_{\nu} \mathcal{R}_{\nu} = (0,1)\times (0,1)$. As an illustration see Fig.(\ref{SupFig4}) where the region $\cup_{\nu}\mathcal{R}_{\nu}$ is plotted when $N$ take values $2$, $5$, $10$, and $100$.

\end{widetext}