\section{Resonator Design}

\begin{figure*}
	\centering
	\includegraphics[width=0.97\linewidth]{figures/device_overview_v5}
	\caption{\textbf{(a)} Cross-section of one period of the FinFET resonator , modeled in COMSOL, with inset highlighting Si fin (dark green) and epitaxial raised source and drain with silicide (blue) \cite{noauthor_comsol_nodate}. Metal and contact layers highlighted in orange and dark gray and dielectric layers left uncolored (light gray). The target mode is defined primarily along the x-direction via the 180 degree phase constraint due to differential drive and sense. \textbf{(b)} Top down view of device shows the central sense transistors with five pairs of capacitors on either side allowing for differential drive, capped on either end by termination (ten gates in the actual devices, abbreviated to three in the figure) and surrounded by a ground ring. First level metal M1 is used for routing within the structure while higher level metals in Mx and Cx (not shown for clarity in this diagram) rest above this routing in either continuous plate or PnC form. The Kx layer represents higher level metal layers used for device to pad routing. \textbf{(c)} Resonator electrical test scheme with DC biasing provided by standalone SMUs, routed through internal PNA bias tees.  }
	\label{fig:deviceoverview}
\end{figure*}

The designed resonator structure, shown in Fig. \ref{fig:deviceoverview}, is fully integrated in the GlobalFoundries 14LPP process and occupies a space of 9 µm x 15 µm \cite{globalfoundries_12lp_2018}. The device is a four port network, with output of the resonator taken as the differential transconductance (Eq. \ref{eq:gmdd}), which is a function of three main items: the electrical to acoustic transduction efficiency in the drive transistors, the acoustic energy confinement and amplification in the resonance cavity, and the acoustic to drain current conversion in the output transistors. The first of these items, drive transduction, is electrostatic due to the 14 nm fabrication process used, although the potential for future piezoelectric drive using CMOS compatible thin film ferroelectrics exists \cite{trentzsch_28nm_2016,he2020tunable,he2019switchable}. The third of these, acoustic-to-drain current conversion, is also largely a function of the process and materials used in device fabrication. While this conversion will be examined utilizing a simplified model, the initial device design is focused on maximizing acoustic energy confinement by varying gate length, termination scheme, and back-end-of-line (BEOL) confinement conditions.

\begin{equation}\label{eq:gmdd}
	g_{mdd} = \frac{i_d}{v_d} = Y_{dd21}-Y_{dd12}
\end{equation}

\subsection{BEOL Phononic Confinement}

Fig. \ref{fig:deviceoverview}(c) depicts a 3D unit cell of the resonator with a 14 nm wide FinFET at the front-end-of-line (FEOL) and periodically arranged metal blocks at the BEOL, categorized into Mx and Cx layers as shown based on thickness and minimum lateral feature sizes. The target resonant mode for the unreleased resonator is along the length direction of the fin transistors. To obtain lower insertion loss and higher quality factor, the radiated acoustic energy needs to be concentrated at the MOS capacitor drive transducers and the sensing transistors. Compared to more traditional suspended resonators, this represents a challenge in the unreleased resonant body transistor design in the commercial CMOS processes.

Integrating a phononic crystal (PnC) at the BEOL of the standard CMOS process and leveraging the mechanical bandgaps of such structures has proven to be an effective tool to confine the energy and boost resonance quality factor \cite{he2020tunable,bahr_theory_2015}. As shown by M1 in \ref{fig:deviceoverview}(b), these periodic metals are uniform along the length (y direction) of the transistor gates (orthogonal to the fin direction) within the resonant cavity.

In order to study and optimize the bandgap induced by the phononic crystals, several 2D simulations were performed in COMSOL Multiphysics. Specifically, by mapping the real coordinates to reciprocal coordinates and searching for the eigenmodes along the first irreducible Brillouin zone in the 2D reciprocal lattice, a 130 nm by 60 nm PnC (such as can be designed within the constraints of the Mx layers in the 14 nm process) yields a calculated bandgap from 8 GHz to 12 GHz (shown in Fig. \ref{fig:DR}(b)).
To further confine the elastic energy, a separate PnC designed under Cx metal layer design rules is also considered for a subset of devices. The corresponding dispersion relationship is shown in Fig. \ref{fig:DR}(c). A mechanical bandgap is then obtained between 7.6 GHz and 10.7 GHz which overlaps with the bandgap induced by the Mx metal PnC.
Together, these PnC designs indicate elastic waves with a frequency of 8 GHz to 12 GHz can be prohibited from propagating into the BEOL. If coupled with a waveguided mode with differential drive that exhibits a large momentum ($k_x = \pi/a$) at a frequency below the sound cone for the bulk silicon substrate, these modes can be confined laterally towards the driving/sensing transistors, with modes from 8 to 10.7 GHz showing the greatest confinement \cite{bahr_theory_2015}. For an in-depth analysis on the acoustic modes in this device, readers are referred to \cite{rawat_analysis_2021}.

\subsection{Resonant Cavity Termination}
 While the PnCs in the main resonant cavity enhance energy confinement and wave reflection from the BEOL, there also exists wave scattering at the two ends of the resonator. Reducing this scattering further increases the horizontal confinement, reduces the insertion loss, and hence increases the Q. The adiabatic theorem in photonic waveguide design is leveraged here, which states that when the wave is propagated down the waveguide, scattering will vanish because of the limitation of sufficiently slow perturbation \cite{meade1995photonic}. Therefore, the key to reduce this scattering is to introduce different units as terminations at the two ends for slower transitions. The unit cell for all the PnC structures in the major termination section throughout the paper is kept as 150 nm x 60 nm.

 While all measured devices have ten terminating gates at each end of the cavity, two types of cavity termination spacings are explored to provide acoustic confinement in the lateral direction. These include abrupt and gradual termination schemes. The abrupt termination refers to the scenario where a terminating periodic array with constant pitch shifted from that of the resonance cavity array is immediately adjacent to the main resonance region. In this case, the gate length of the termination immediately transitions to the $L_{term}$ value given in Fig. \ref{fig:rfoverview}(a) and stays  constant until the end of the device. As such, the guided waves are abruptly transitioned from the main cavity towards the ends. In the gradual termination scenario, the main resonant cavity periodicity adiabatically transitions from $L_{gate}$ to the terminating array periodicity $L_{term}$, resulting in a termination gate length gradually transitioning from 80 to 140 nm. This serves to reduce scattering from the two ends and boost the Q \cite{bahr2015theory}. These two schemes are depicted in Fig. \ref{fig:rfoverview}(c).

\subsection{FinFET Sensing}

To understand the conversion from stress to drain current in the sense transistors, modulation of three separate FET properties is examined: oxide capacitance, channel mobility due to piezoresistivity, and silicon bandgap.
For ease of visualization, an initial analysis is performed with the source-referenced simplified strong-inversion model with $\alpha=1$ (the traditional "square-law model", or SPICE level 1) \footnote{This model has only a handful of parameters to describe device operation, as opposed to the more advanced BSIM models used commonly in industry, which include many empirically fit parameters to give the best device models for commercial applications \cite{noauthor_bsim_nodate}. This simplified model is used with the understanding that the derivation, while not capturing a variety of short channel effects that influence modern transistor behaviour, will provide a more intuitive understanding into the basics of active FET sensing.} \cite{tsividis_operation_2010}:

\begin{equation}
	I_d = \mu C_{ox} \frac{W}{L}\Big((V_{GS}-V_T)V_{DS} - \frac{V_{DS}^2}{2}\Big)
\end{equation}
where $I_d$ is the drain current, $\mu$ is the carrier mobility, $C_{ox}$ is the gate capacitance, $V_{DS}$ and $V_{GS}$ are the drain-to-source and gate-to-source bias, respectively, $W$ and $L$ are the channel width and length, and $V_T$ is the threshold voltage required for significant channel conduction.

Capacitance modulation in these devices is the same mechanism used in electrostatic drive of these devices, except that rather than sensing the gate current directly, the capacitance modulation is being amplified through transistor action into a change in drain current. The piezoresistive effect, whereby the semiconductor carrier mobility is modulated through stress in the transistor channel, likewise appears as a multiplier on the transistor drain current. Semiconductor bandgap modulation, on the other hand, largely effects the threshold voltage of the transistor. For a complete derivation of these effects, the reader is directed to the derivation in Appendix \ref{sec:derivation}.
Modified to include the effects of stress on $\mu$, $C_{ox}'$, and $V_T$ (modeled as $\mu(\sigma) = \mu_0+\Delta \mu$), the original expression for current in the linear region becomes the following:

\begin{multline} \label{eq:id_stressfet}
	I_d(\sigma)	= (\mu_0 C_{ox0}' +\mu_0 \Delta C_{ox}' + C_{ox0}'\Delta \mu + \Delta \mu\Delta C_{ox}')  *\frac{W}{L}\Big((V_{GS}-V_{T0})V_{DS} - \frac{V_{DS}^2}{2}\Big) \\
	-\Delta V_T V_{DS} \frac{W}{L} (\mu_0 C_{ox0}' +\mu_0 \Delta C_{ox}' + C_{ox0}'\Delta \mu + \Delta \mu\Delta C_{ox}')
\end{multline}
where
\begin{multline} \label{eq:deltaid_stressfet}
	\Delta I_d(\sigma)	= (\mu_0 \Delta C_{ox}' + C_{ox0}'\Delta \mu + \Delta \mu\Delta C_{ox}') *\frac{W}{L}\Big((V_{GS}-V_{T0})V_{DS} - \frac{V_{DS}^2}{2}\Big) \\
	-\Delta V_T V_{DS} \frac{W}{L} (\mu_0 C_{ox0}' +\mu_0 \Delta C_{ox}' + C_{ox0}'\Delta \mu + \Delta \mu\Delta C_{ox}')
\end{multline}

Likewise, in saturation

\begin{equation}
	I_{d0} = \mu C_{ox} \frac{W}{2L}\Big((V_{GS}-V_T)^2\big(1+\lambda[V_{DS}-(V_{GS}-V_{T0})]\big)\Big)
\end{equation}
where $\lambda$ represents the channel length modulation coefficient (0.05 to 0.25 typical) and

\begin{multline}
	\Delta I_d
	= (\mu_0 \Delta C_{ox}' + C_{ox0}'\Delta \mu + \Delta \mu\Delta C_{ox}')
	*\frac{W}{2L}\Big((V_{GS}-V_T)^2(1+\lambda(V_{DS} - V_{GS} + V_{T}))\Big) \\
	\hspace{13mm}- \frac{W}{2L} (\mu_0 C_{ox0} +\mu_0 \Delta C_{ox} + C_{ox0}\Delta \mu + \Delta \mu\Delta C_{ox}) \Bigg[ \hfill \\
	\hspace{27mm}\Delta V_T \Big[\lambda(V_{GS}-V_{T0})^2 - 2(V_{GS}-V_{T0})
	*\Big(1+\lambda\big[V_{DS}-(V_{GS}-V_{T0})\big]\Big)\Big] \hfill \\
	\hspace{23mm}+ \Delta V_T^2 \Big[\Big(1+\lambda\big[V_{DS}-(V_{GS}-V_{T0})\big]\Big) - 2\lambda(V_{GS}-V_{T0})\Big] \hfill \\
	\hspace{23mm}+ \Delta V_T^3 [\lambda]  \Bigg] \hfill
\end{multline}

For modern sub-micron devices, carriers do not exhibit constant mobility but rather undergo velocity saturation at high fields. A simplified way to model this analytically is through a piece-wise equation
\begin{equation}
	v_d = \begin{cases}
		\frac{\mu |E_x|}{1+|E_x|/E_{sat}} \quad,\quad |E_x| < E_{sat} \\
		v_{sat}\hfill,\quad |E_x| > E_{sat}
	\end{cases}
\end{equation}
which, solved at $v_d = v_{sat}$, gives
\begin{equation} \label{eq:esat}
	E_{sat} = \frac{2v_{sat}}{\mu}
\end{equation}
and leads to a modified expression for drain current
\begin{equation}
	I_d = \begin{cases}
		\mu C_{ox} \cfrac{W}{L}\cfrac{1}{1+\cfrac{V_{DS}}{E_{sat}L}}\Big((V_{GS}-V_T)V_{DS} - \frac{V_{DS}^2}{2}\Big) \quad,\quad V_{DS} < V_{DSsat} \\ \\
		C_{ox} W v_{sat}\cfrac{(V_{GS}-V_T)^2}{(V_{GS}-V_T)+E_{sat}(L-\Delta L)} \hfill,\quad V_{DS} > V_{DSsat}
	\end{cases}
\end{equation}
where
\begin{equation}
	V_{DSsat} = \frac{(V_{GS}-V_T)E_{sat}L}{(V_{GS}-V_T)+E_{sat}L}.
\end{equation}

An important nuance for acoustic resonators is that both the saturation velocity and the change in mobility (piezoresistive effect) are related to the change in carrier effective mass caused by acoustic deformation of the crystal lattice - saturation velocity by $m^{*-1/2}$ and mobility by a factor between $m^{*-1/2}$ and $m^{*-5/2}$, depending on doping (for a detailed analysis, see Appendix \ref{sec:mobility}). This relation means that, while DC current may saturate, the AC piezoresistive effect can still change drain current to an extent, even in the saturation region. For this reason, a compound model is utilized that incorporates velocity saturation in the calculation of DC drain current but allows mobility to vary freely. Under this approach, the contribution to RF output from both electrostatic forces and mobility modulation should scale linearly with drain current. The contribution from any threshold voltage shift is seen to be related to drain bias only in the linear regime and is scaled by the impact of mobility and capacitance changes. In the saturation regime, on the other hand, the threshold voltage shift is primarily a function of gate bias and is modified by drain bias only through channel length modulation. These expressions for $\Delta I_d$ represent one half of the device transfer function $g_{mdd}$, discussed in Eq. \ref{eq:gmdd}, with the other half being given by the voltage-to-strain conversion and subsequent strain confinement/resonance. While direct measurement of strain in the solid-state resonator in the front end of line poses a challenge, all three of these aspects are subsequently investigated using measured electrical results.

\section{Results and Discussion}

DC results, shown as an example for the 80 nm gate length device A1 in Appendix \ref{sec:sup_figs}, indicate that the sense transistors in all devices were working as intended. The large current values of approximately 2 mA when gate and drain are biased at $V_{dd}$ are due to the parallel DC connection of the two sense transistors and large effective width of the individual transistors, where 40 fins span the resonance cavity (in y-direction) and are electrically connected in parallel.

\subsection{Analysis of Design Variations}

\begin{figure*}
	\centering
	\includegraphics[width=1\linewidth]{figures/sample_overview_v5}
	\caption{\textbf{(a)} Depiction of the termination schemes investigated for use at the ends of the resonant cavity. \textbf{(b)} Phononic crystal structure, with various investigated pitches. \textbf{(c)} Device design variations and \textbf{(d)} their experimentally measured differential transconductance at $V_D=0.6$ V, $V_G=0.8$ V, $V_{DR}=0.8$ V. \textbf{(e)} Target mode, simulated in COMSOL, with stress concentrated in finFET channel and confined by the BEOL metal layers.}
	\label{fig:rfoverview}
\end{figure*}

While DC results were similar for all devices (outside expected differences due to gate length), resonance mode confinement was highly dependent on the design variations studied. A summary of all measured devices as well as a comparison of their differential transconductance spectra from 9.5 to 13.5 GHz are shown in Fig. \ref{fig:rfoverview}. From these data, it is evident that phononic crystal confinement in the Mx layers results in significantly improved quality factor and maximum transconductance. Also visible is the role the type of termination has in defining the center frequency of nearby spurious modes in the device.

\begin{table*}[]
	\renewcommand{\arraystretch}{1.3}
	\caption{ANOVA Results for Peak Fits}
	\label{tab:anova_results}
	\centering
	\resizebox{\textwidth}{!}{%
	\begin{tabular}{|c|c|c||c|c|c|c|c|c||c|}
		\hline
		\multicolumn{2}{|c|}{\scriptsize{95\% Confidence Intervals} } & \multirow{2}{*}{n} & \multicolumn{2}{c|}{Frequency [Hz]} & \multicolumn{2}{c|}{Height [nS]} & \multicolumn{2}{c||}{Q} & \multirow{2}{*}{Samples} \\
		\cline{4-9}
		\multicolumn{2}{|c|}{\scriptsize{Significant results in bold.}} &  & Mean & Std. Error & Mean & Std. Error & Mean & Std. Error &  \\
		\hhline{|==|=||=|=|=|=|=|=||=|}
		\rule[-1ex]{0pt}{2.5ex} \multirow{3}{*}{PnC Confinement} & Mx, Cx & 413 & 11.27 G & 52.8 M & \textbf{451.1} & \textbf{30.6} & \textbf{44.32} & \textbf{1.68} & \multirow{3}{*}{\shortstack{A1-A4, B0-B1\\ (L$_{gate}$=80 nm)}} \\
		\cline{2-9}
		\rule[-1ex]{0pt}{2.5ex}  & Mx Only & 368 & 11.22 G & 55.9 M & \textbf{379.4} & \textbf{32.5} & \textbf{37.55} & \textbf{1.78} &  \\
		\cline{2-9}
		\rule[-1ex]{0pt}{2.5ex}  & Plates & 276 & 11.26 G & 64.6 M & \textbf{203.6} & \textbf{37.5} & \textbf{19.88} & \textbf{2.05} &  \\
		\hline
		\rule[-1ex]{0pt}{2.5ex} \multirow{2}{*}{Termination} & Abrupt & 1289 & \textbf{11.16 G} & \textbf{33.4 M} & \textbf{301.6} & \textbf{13.9} & \textbf{24.21} & \textbf{0.795} & \multirow{2}{*}{All} \\
		\cline{2-9}
		\rule[-1ex]{0pt}{2.5ex}  & Gradual & 1145 & \textbf{11.29 G} & \textbf{35.4 M} & \textbf{352.4} & \textbf{14.7} & \textbf{ 31.71} & \textbf{0.843} &  \\
		\hline
		\rule[-1ex]{0pt}{2.5ex} \multirow{2}{*}{L$_{gate}$} & 80 nm & 30 & \textbf{11.66 G} & \textbf{23.0 M} & 1,866 & 153 & 84.8 & 4.40 & \multirow{2}{*}{\shortstack{A3-A6 \\ Q$>$25, Height$>$500 nS}} \\
		\cline{2-9}
		\rule[-1ex]{0pt}{2.5ex}  & 106 nm & 30 & \textbf{10.71 G} & \textbf{23.0 M} & 1,651 & 153 & 87.1 & 4.40 &  \\
		\hline
	\end{tabular}}
\end{table*}

These differences are quantified in Table \ref{tab:anova_results}, where the center frequency, quality factor, and height of resonant peaks in $g_{mdd}$ are compared for the three design splits after elimination of high-error peak fit ($\approx$5\% of data overall). For termination, this process is straightforward as there exists an abrupt and gradually terminated copy of every device. The resulting analysis indicates an increased transconductance in gradual transition devices (352 nS) vs. abruptly terminated ones (302 nS), with a corresponding increase in average quality factor (32 vs. 24). In addition to this enhancement, average resonance center frequency was shown to decrease slightly from 11.29 to 11.16 GHz, driven by the downward shift in the left-most peak of the high-$Q$ cluster of modes seen in PnC-confined devices.

To investigate phononic crystal confinement, analysis was limited to designs with 80 nm gate length (A1-A4 and B0-B1) and the confinement treated as a three-level categorical factor. This analysis showed no significant change in center frequency of the modes, but it did show significantly improved transconductance (204, 379, 451 nS) and quality factor (20, 38, 44) with increasing layers of phononic crystal confinement. These trends are highlighted in Fig. \ref{fig:peakfitresults} along with an example of the peak fitting results used for parameter extraction.

\begin{figure*}
	\centering
	\includegraphics[width=0.9\linewidth]{figures/peakfit_results_v2}
	\caption{ \textbf{(a)}: Output of peak fit algorithm for one of the 392 individual measurements evaluated. Peak fits with error above 0.025 were filtered (120 out of 2554 peaks) and then design splits were evaluated by comparing means and performing ANOVA between appropriate samples. The most significant process variation, BEOL confinement, is shown in \textbf{(b)}, with complete results given in Table \ref{tab:anova_results}. In addition to increasing average peak height and $Q$ across a range of drive and drain biases, it is also evident that an Mx-PnC greatly increases the maximum peak height, from under 1 µS to 4.5 µS for the highest-performing resonator. \textbf{(c)} shows Q and amplitude of fit resonance peaks, filtered by a minimum $Q$ of 25 and height of 750 nS, with each dot indicating a peak of the correspondingly colored device under a single bias condition. Transconductance variations at the same frequency show the dependence of the resonance on electrical biases. Also evident is the change in center frequency due to gate length, and splitting from two to three distinct modes with gradual termination, as first shown in Fig. \ref{fig:rfoverview}. }
	\label{fig:peakfitresults}
\end{figure*}

Compared to the other factors studied, L$_{gate}$ is unique in that it is the primary design variable for selecting center frequency due to the mode being defined by source/drain spacing. This means that the different values studied dramatically change the cavity shape and thus all peaks detected in the studied frequency span of 9.5 to 13.5 GHz cannot be used in the analysis, as they are not present in this range in all devices. To overcome this, a much smaller comparison of A3 and A4 (L$_{gate}$=80 nm) to A5 and A6 (L$_{gate}$=106 nm) was performed using only the prominent modes with transconductance above 500 nS and $Q$ above 25 that were present towards the center of the spectra. The downside of this approach is that gate length is confounded with a change in phononic crystal design that is necessary for the different device geometries due to design rule challenges. Despite this, the fact that confinement method showed no significant change in center frequency for 80 nm devices points towards this shift being dominated by gate length as would be expected.

As the device with the highest performance, resonator A2 was selected for further analysis, with the maximum transconductance of the $\approx$11.75 GHz mode examined as a function of device biasing. This was done at a drive bias of $V_{dd}$, as $g_{mdd}$ was seen to increase linearly with the increase in MOS capacitance at higher drive bias (see Appendix \ref{sec:sup_figs}). The dependence of $g_{mdd}$ on sense transistor biasing, shown in Fig. \ref{fig:exp_vs_sim}, highlights that device transconductance, as expected, is modulated in a manner consistent with piezoresistive and bandgap modulation. While transistor parameters are of the correct order of magnitude in the reported data fit, they vary slightly from real-world values due to the simplified model used. For example, the square law model does not accurately capture subthreshold (accumulation and weak inversion) behaviour, and thus the fit favors a slightly lower Vt to more closely fit bias currents near zero. This, in turn, leads to an overestimation of current at higher biases, which is compensated by a lower than expected value for mobility and a slightly larger equivalent oxide thickness. On the AC side, the fit is hindered by the presence of a noise floor in the s-parameter measurements that may obscure lower peak levels. This, along with the poor near-threshold characteristics of the square law model, contributes to the fitting error at lower values of $g_{mdd}$. Moving away from the square law model to a more advanced model such as that derived by Taur et al. that captures all regions of operation may decrease this error \cite{taur_continuous_2004}. Likewise, adopting band-structure simulations of silicon as a result of acoustic deformation could provide a unified source of mobility and threshold voltage changes in the transistor, at the expense of increased computation requirements.

For NMOS devices oriented in the <110> direction, the negative piezoresistive coefficients ensure that the contribution of the threshold voltage shift and mobility is additive. If fRBT devices are implemented with PMOS devices or in processes with different fin orientations, care should be taken to ensure that the dominant stress-dependent effects do not compete against each other. This may not be as much of an issue for long channel RBT devices, as piezoresistive modulation dominates for long channel devices that are not supply voltage constrained.

\begin{figure*}
	\centering
	\includegraphics[width=1\linewidth]{figures/exp_vs_sim.png}
	\caption{\textbf{(a)} Simulated contributions to $\Delta I_D$ from bandgap modulation (threshold voltage), piezoresistive modulation (mobility), electrostatic modulation (oxide capacitance), and cross-term. For short channel devices that are supply voltage constrained, bandgap and piezoresistive modulation are of similar magnitudes and dominate over electrostatic terms.  \textbf{(b)} These two terms, along with the overall drain current modulation, are plotted against DC bias current of the sense transistor (X-axis) for simulated (top) and experimental (bottom) devices. Gate bias is shown by marker color (blue to yellow) and ranges from 0.4 to 0.8 by 0.1 V. Drain bias is represented with marker size, moving from small to large as bias increases from 0.1 to 0.8 by 0.1 V. The large cluster of points with equal gate biases and similar drain currents indicate saturation mode operation for those bias conditions with channel length modulation. The mobility term is linear with drain current, whereas the threshold voltage term due to bandgap modulation has both a flat \textbf{(i)} region in the linear regime with no $V_G$ dependence for constant $V_D$ and a slight $V_D$ dependence \textbf{(ii)} for constant $V_G$ in the saturation regime. Experimental data shows roughly linear dependence of drain current on gate voltage (a hallmark of velocity saturation \cite{tsividis_operation_2010} vs. the $V_G^2$ dependence of the long channel model). Some deviation between the model fit and experimental results can be seen at high drain bias and low gate bias, which may be indicative of mobility degradation due to vertical fields, which is not modeled currently. Otherwise, experimental results show a similar shape to that of the simulated device, indicating the dominance of both piezoresistive and bandgap modulation in the experimental resonator. \textbf{(c)} Model fit to experimental results ($\mu=79$ cm$^2$/Vs, equivalent oxide thickness=1.6 nm, $V_t$=0.206 V, $v_{sat}$=1.83*10$^7$ cm/s, $\lambda$=0.47, $W$ = 3 um, total stress=1 MPa). A basin-hopping algorithm is used to find a global minimum fit for the DC transistor parameters, which are then fixed for the AC fit, achieved by varying stress.}
	\label{fig:exp_vs_sim}
\end{figure*}

The highest magnitude peak fit occurred in device A2 (L$_{gate}$ of 80 nm, gradually terminated, dual PnC) at 11.73 GHz with a differential transconductance of 4.49 µS and a $Q$ of 69.8, for an $f*Q$ product of 8.19*10$^{11}$ at drain, gate, and drive biases of 0.6, 0.8, and 0.8 V, respectively. While previous unreleased resonators in the 32 SOI process have achieved an $f*Q$ product of $3.8*10^{13}$ with phononic confinement at 3.3 GHz, these 14 nm devices exhibit 50x higher transconductance and show an almost 3x improvement in Q over a similar frequency measured mode in 32 SOI devices \cite{bahr_optimization_2016, marathe_resonant_2014}. Compared to some contemporary resonators in the X and adjacent K$_\textrm{u}$ (12-18 GHz) bands in Table \ref{tab:lit_comparison}, the unreleased devices in this work achieve moderate $f*Q$ in a commercially available process in an area smaller than comparable released devices. This, combined with the fact that the devices are integrated in the front-end-of-line of a commercial CMOS advanced technology node, opens the door for tight logic integration with minimum routing (and associated parasitics). Additionally, if process design rule challenges can be overcome, this performance can be achieved across a band of frequencies by varying gate length and phononic crystal design.

\begin{table*}
	\renewcommand{\arraystretch}{1.3}
	\centering
	\caption{Comparison with Contemporaneously Reported X and K$_\textrm{u}$ Band Resonators}
	\label{tab:lit_comparison}
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
		Ref & $f_0$ [GHz] & $Q$ & $f_0 * Q$ & Technology & Released & Area [um$^2$] \\
		\hline
		\textbf{This Work} & 11.7 & 69.8 & 8.2*10$^{11}$ & 14LPP & No & 135 \\
		\hline
		\cite{marathe_resonant_2014} & 11.5 & 24 & 2.8*10$^{11}$ & 32SOI  & No & 15 \\
		\hline
		\cite{chen_high-q_2019} & 8.8 & 750 & 6.6*10$^{12}$ & Custom AlN & Yes & 6,500 \\
		\hline
		\cite{assylbekova_11_2020} & 11.1 & 615 & 6.8*10$^{12}$ & Custom AlN & Yes & 675 \\
		\hline
		\cite{yang_1060-ghz_2020} & 13.0 & 282 & 3.7*10$^{12}$ & Custom LiNbO3 & Yes & 2,000 \\
		\hline
		\cite{abouyoussef_quad_2018} & 14.5 & 400 & 5.8*10$^{12}$ & PCB Microstrip  & - & 27,000,000 \\
		\hline
	\end{tabular}
\end{table*}

\section{Conclusion}

This work provides an in-depth experimental investigation into the influence of several design parameters on the performance of unreleased Resonant Fin Transistors in standard 14nm FinFET technology. Phononic crystal confinement, shown experimentally here to be the most important variation in improving unreleased resonator performance (2.2x improvement on average), provides significant opportunity for integrating acoustic devices into standard CMOS platforms. The highest performing resonator, a gradually terminated 80 nm gate length device with both Mx and Cx phononic crystal confinement, shows a transconductance amplitude and $Q$ of 4.49 µS and 69.8 respectively, with an $f*Q$ product of 8.19*10$^{11}$. At the same time, an analytical model has been demonstrated that predicts resonator performance across sense transistor biasing conditions, utilizing a combined velocity saturated current at DC and unsaturated mobility variation at AC, where acoustic perturbation of carrier effective mass modulates saturation velocity and mobility simultaneously. This highlights the importance of understanding the assumptions that go into developing advanced electrical models for short channel transistors, which may not accurately depict effects present in active electromechanical devices. While creating these phononic confinement schemes in a standard process with strict design rules poses a unique challenge, we have established a methodical way to predict and design resonator performance across multiple CMOS platforms. This technology combined with parallel advancements in CMOS-compatible ferroelectric and piezoelectric thin films, presents the opportunity to develop compact, low-cost, CMOS-integrated, and electrically-controllable resonators with no additional packaging required, opening the door for acoustic processing alongside traditional electronics in advanced CMOS nodes.


\section{Introduction}

The existence of gravitational waves (GWs), considered as ripples in space-time fabric,
was predicted in Einstein's theory of general relativity in 1916.
It was only until recently that GWs emitted by a binary black hole merger were detected by the LIGO experiment~\cite{ligo}.
Since the first detection, many GW events from compact binary coalescences have been observed in Advanced LIGO and Advanced Virgo O2 and O3 observing runs~\cite{ligo2019}. This opens up a new era of GW astronomy and cosmology. Upcoming and future GW experiments such as KAGRA, GEO600, Einstein Telescope, Cosmic Explorer, LISA, DECIGO, Taiji, and TianQin, will bring us a precision science in GW observation~\cite{ligo2050}.

Stochastic gravitational wave background (SGWB) is one of the main goals in GW experiments. There have been many studies on potential astrophysical and cosmological sources for SGWB such as distant compact binary coalescences, early-time phase transitions, cosmic string or defect networks, second-order primordial scalar perturbations, and inflationary GWs~\cite{romano}. GWs are very weakly interacting, so they decouple from matter at the time of production and then travel to us almost without being scattered. At the present, they remain as a SGWB that carries initial conditions of the production processes in the early Universe.

In general, the SGWB can be anisotropic and polarized. For examples, helical GWs can be produced in axion inflation models, leading to a net circular polarization~\cite{alexander,satoh,sorbo,crowder2013}. Linear polarization can be generated through diffusion by compact astrophysical objects, with an amount suppressed by a factor of at least $10^{-4}$ with respect to the intensity anisotropies; however, it can be enhanced if dark matter is dominated with sub-solar-mass primordial black holes~\cite{cusin}. Furthermore, the directionality dependence of the SGWB have been recently explored~\cite{bartolo,pitrou}.

The method adopted in current GW interferometry such as the LIGO-Virgo experiment for detecting the intensity and polarization anisotropies or the Stokes parameters of the SGWB is to correlate the responses of a pair of detectors to the GW strain amplitude. The correlation allows us to filter out detector noises and increase the signal-to-noise ratio~\cite{romano}. The on-going and future GW interferometry experiments have sensitivities for detecting GWs at frequencies ranging from kilohertz to millihertz.

Another method to search for SGWB anisotropies is through their gravitational effects on the arrival times of radio pulses from millisecond pulsars~\cite{romano}. In the pulsar-timing observation, radio pulses from an array of roughly 100 Galactic millisecond pulsars are being monitored with ground-based radio telescopes. The fluctuations in the redshifts of the radio pulses due to the presence of GWs would allow us to detect longer-wavelength GWs at nanohertz frequencies. Current pulsar timing array (PTA) experiments include EPTA~\cite{EPTA},  NANOGrav~\cite{NANO}, and PPTA~\cite{PPTA}. The future SKA radio telescope will have projected sensitivity three to four orders of magnitude better than the current PTAs~\cite{SKA}.
Recently, the NANOGrav Collaboration~\cite{nanograv} has found strong evidence of a stochastic common-spectrum process across 45 millisecond pulsars, alluding to a SGWB with a characteristic strain of $h_c=1.92\times10^{-15}$
at a reference frequency of $f_{\rm yr}=1\,{\rm yr^{-1}}\simeq 31.8\,{\rm nHz}$.
However, they have not found statistically significant evidence that this process has quadrupolar spatial correlations, as expected for the properties in the presence of a SGWB.

The pulsar-timing observation of the SGWB has been studied extensively~\cite{romano}. The fluctuation in the redshift of the radio pulse from a pulsar in the presence of GWs is due to the Shapiro time delay given by the time-varying gravitational potential difference from the pulsar to the Earth. By ignoring the contribution at the pulsar site (the so-called pulsar term), the overlap reduction functions (ORFs) for the SGWB intensity anisotropy have been derived using the two-point correlation method~\cite{anholm,mingar13}. Later, it has been extended to including the SGWB circular polarization anisotropy~\cite{kato16}. Furthermore, the corrections of the ORFs for the SGWB intensity anisotropy by the pulsar term have been discussed~\cite{mingar14}. In this paper, we develop a numerical scheme to compute the ORFs for the SGWB intensity and polarization anisotropies, following the spherical harmonic expansion method adopted in Ref.~\cite{chu21} though in the context of GW interferometry.  Instead of using the Shapiro time delay, we follow Ref.~\cite{ng21} to begin with the Sachs-Wolfe line-of-sight integral and then calculate the two-point correlation function. This scheme allows us to compute the ORFs for high multipole moments order by order and simultaneously incorporates the effects of the pulsar term.

\section{Polarized SGWB}

In the Minkowskian vacuum, the metric perturbation $h_{ij}$ in the transverse traceless gauge depicts GWs propagating at the speed of light $c=\omega/k$. At a given spacetime point $(t,\vec{x})$, it can be expanded in terms of its Fourier modes:
\begin{equation}
\label{eq:planwave}
h_{ij}(t,\vec{x}) = \sum_{A}\int_{-\infty}^{\infty} \d f \int_{S^2} \d\hat{k} \;
    h_A(f,\hat{k}) \mathbf{e}^A_{ij}(\hat{k})
    e^{-2 \pi i f (t - \hat{k}\cdot \vec{x}/c)}\,,
\end{equation}
where $A$ stands for the polarization or the helicity of GWs described by the corresponding basis tensors $\mathbf{e}^A_{ij}(\hat{k})$, which are transverse to the direction of the wave propagation denoted by $\hat{k}$. Since $h_{ij}$ is real, its Fourier components are not fully independent with each other. For our application, we require those Fourier components with negative frequencies to be $h_A(-f,\hat{k})=h_A^*(f,\hat{k})$ for all $f\ge 0$. The GWs are considered as stochastic as long as $h_{ij}$ are random fields thus characterized by their ensemble averages. Besides, assuming the probability distribution of the random amplitude $h_{ij}$ be Gaussian, then only the two-point correlation function $\langle h_{ij}(t,\vec{x}_1) h_{ij}(t,\vec{x}_2) \rangle$ is needed to describe its statistical behavior. Furthermore, if the waves are homogeneous, i.e., having translational symmetry, the ensemble average can be evaluated by doing spatial averages. As a result, the two-point correlation function of the Fourier modes should have the following form
\begin{equation}
\label{eq:paa}
\langle h_{A}(f,\hat{k}) h^*_{A'}(f',\hat{k}') \rangle
= \delta(f-f') \delta(\hat{k}-\hat{k}')P_{AA'}(f,\hat{k}) \,,
\end{equation}
where the $\delta(f-f')$ arises from the delta function of the magnitude of their 3-momenta $\delta(\vec{k}-\vec{k}')$ and the assumption made in Eq.~(\ref{eq:planwave}) that these waves satisfy the equation of motion in vacuum. Also, the presence of $\delta(f-f')$ implies that the signal is stationary. This is a fairly good approximation during an observing period for a typical experiment.

For GWs coming from the sky direction $-\hat{k}$ with wave vector $\vec{k}$, it is customary to write the polarization basis tensors in terms of the basis vectors in the spherical coordinates:
\begin{align}
    \mathbf{e}^{+}(\hat{k}) &= \hat{\mathbf{e}}_\theta \otimes \hat{\mathbf{e}}_\theta
                               - \hat{\mathbf{e}}_\phi \otimes \hat{\mathbf{e}}_\phi \,, \nonumber \\
    \mathbf{e}^{\times}(\hat{k}) &= \hat{\mathbf{e}}_\theta \otimes \hat{\mathbf{e}}_\phi
                                    + \hat{\mathbf{e}}_\phi \otimes \hat{\mathbf{e}}_\theta  \,,
\end{align}
in which $\hat{\mathbf{e}}_\theta$, $\hat{\mathbf{e}}_\phi$, and $\hat{k}$ form a right-handed orthonormal basis.
Also, we can define the complex circular polarization basis tensors as
\begin{align}
    \mathbf{e}_{R} &= \frac{(\mathbf{e}_{+} + i \mathbf{e}_{\times})}{\sqrt{2}} \,,
   &\mathbf{e}_{L} &= \frac{(\mathbf{e}_{+} - i \mathbf{e}_{\times})}{\sqrt{2}} \,,
\end{align}
where $\mathbf{e}_{R}$ stands for the right-handed GW with a positive helicity while $\mathbf{e}_{L}$ stands for the left-handed GW with a negative helicity. The corresponding amplitudes in Eq.~(\ref{eq:planwave}) in the two different bases are related to each other via
\begin{align}
    h_{R} &= \frac{(h_{+} - i h_{\times})}{\sqrt{2}} \,,
   &h_{L} &= \frac{(h_{+} + i h_{\times})}{\sqrt{2}} \,.
\end{align}

Analogous to the case in electromagnetic waves~\cite{book:BornAndWolf}, the coherency matrix $P_{AA'}$ in Eq.~(\ref{eq:paa}) is related to the Stokes parameters, $I$, $Q$, $U$, and $V$ as

\begin{align}
    I &= \left[ \langle h_R h_R^* \rangle + \langle h_L h_L^* \rangle \right] / 2 \,,\nonumber\\
    Q + iU &=  \langle h_L h_R^*  \rangle \,,\nonumber\\
    Q - iU &=  \langle h_R h_L^*  \rangle \,,\nonumber\\
    V &= \left[ \langle h_R h_R^* \rangle - \langle h_L h_L^* \rangle \right] / 2 \,\label{eq:spIQUV}.
\end{align}
They are functions of the frequency $f$ and the propagation direction $\hat{k}$.
To get some flavor of the meaning of these Stokes parameters, we may take a look at an example for an unpolarized quasimonochromatic GW signal with a constant intensity. It should have a constant $I$ which represents the total intensity regardless of its polarization. We have $V=0$ since the power in the right-handed and the left-handed modes should be identical. Also, because the relative phase between the left-handed and the right-handed modes ($\arg(h_L)-\arg(h_R)$) is random for an unpolarized source, the ensemble average $\langle h_L h_R^* \rangle \sim \langle e^{i(\arg(h_L)-\arg(h_R))} \rangle$ becomes zero, thereby making $Q=U=0$.

\section{Pulsar timing}
\label{ptiming}

In the pulsar-timing observation, radio pulses from an array of roughly 100 Galactic millisecond pulsars are being monitored with ground-based radio telescopes. The redshift fluctuation of a pulsar in the pointing direction
$\hat{e}$ on the sky is given by the Sachs-Wolfe effect~\cite{sachs},
\begin{equation}
z(\hat{e})= - {1\over 2}\int_{\eta_e}^{\eta_r}\d\eta\, \hat{e}^i \hat{e}^j \frac{\partial}{\partial\eta}h_{ij}(\eta, \vec x)\,,
\label{rsfze}
\end{equation}
where the lower (upper) limit of integration in the line-of-sight integral represents the point of emission (reception)
of the radio pulse. The physical distance of the pulsar from the Earth is
\begin{equation}
D=c(\eta_r-\eta_e)\,,
\end{equation}
which is of order $1\,{\rm kpc}$.

The quantity that is actually observed in the pulsar-timing observation is the time residual counted as
\begin{equation}
r(t)=\int_0^t dt'z(t')\,,
\end{equation}
where $t'$ denotes the laboratory time and $t$ is the duration of the observation.
Using the laboratory time $t'$, we rewrite Eq.~(\ref{rsfze}) as
\begin{equation}
z(t',\hat{e})= - {1\over 2}\int_{t'+\eta_e}^{t'+\eta_r}\d\eta\, d^{ij}\frac{\partial}{\partial\eta}h_{ij}(\eta, \vec x)\,,
\end{equation}
where the detector tensor is
\begin{equation}
d^{ij}=\hat{e}^i \hat{e}^j\,.
\end{equation}
Then, replacing $\vec{x}$ by $c(\eta_r-\eta)\hat{e}$, we have
\begin{widetext}
\begin{align}
r(t)&={1\over2}\sum_{A}\int_{-\infty}^{\infty} \d f \int_{S^2} \d\hat{k} \int_{\eta_e}^{\eta_r}\d\eta\,(1-e^{-2\pi i f t}) \,
    h_A(f,\hat{k}) \,d^{ij} \mathbf{e}^A_{ij}(\hat{k})\,
    e^{-2 \pi i f \eta} e^{2\pi i f (\eta_r-\eta) \hat{k} \cdot \hat{e}} \nonumber \\
    &=\frac{1}{4\pi}\sum_{A}\int_{-\infty}^{\infty} \frac{\d f}{f} \int_{S^2} \d\hat{k}\, (1-e^{-2\pi i f t})\, e^{2\pi i f\eta_r}
    h_A(f,\hat{k})\, d^{ij} \mathbf{e}^A_{ij}(\hat{k}) \left[1-e^{2\pi i f D (1+\hat{k} \cdot \hat{e})/c}\right]\frac{i}{1+\hat{k} \cdot \hat{e}}\,.
\end{align}
\end{widetext}
This is the result for the time residual due to the Shapiro time delay, where the first term in the square brackets is called the Earth term and the second is the pulsar term. The last term has a pole at $\hat{k} \cdot \hat{e}=-1$. For Galactic pulsars in a pulsar-timing experiment with a typical duration of observation, we have $fD\gg 1$, so the pulsar term is a highly oscillatory function in the $\hat{k}$-integral. If the separation angle between any pair of pulsars is large enough, we can safely omit the pulsar term to simplify the $\hat{k}$-integral for obtaining the two-point correlation function~\cite{anholm,mingar13,kato16}.

In the present work, we retain the Sachs-Wolfe line-of-sight integral. This allows us to avoid the pole term and do the spherical harmonic expansion. In addition, the effect of the pulsar term is now put back in the time integral. Then, using the spherical wave expansion~(\ref{eq:swexpansion}) for the phane wave~(\ref{eq:planwave}), we obtain
\begin{widetext}
\begin{equation}
r(t)=2\pi\sum_{A}\int_{-\infty}^{\infty} \d f \int_{S^2} \d\hat{k} \int_{\eta_e}^{\eta_r}\d\eta\,(1-e^{-2\pi i f t})
    h_A(f,\hat{k}) d^{ij} \mathbf{e}^A_{ij}(\hat{k})
    e^{-2 \pi i f \eta} \sum_{LM} i^L j_L[2\pi f (\eta_r-\eta)] Y_{LM}^*(\hat{k}) Y_{LM}(\hat{e})\,.
\end{equation}
\end{widetext}

\section{Residual Correlation}

The time-residual correlation between a pair of pulsars $a$ and $b$ is constructed as
\begin{equation}
\langle r(t_a)r(t_b) \rangle = \int_0^{t_a}dt' \int_0^{t_b}dt{''} \langle z(t',\hat{e}_a) z(t{''},\hat{e}_b)\rangle.
\end{equation}
Using Eq.~(\ref{eq:paa}) and defining $x=2\pi f (\eta_r-\eta)$, this becomes
\begin{widetext}
\begin{align}
\label{eq:xi2gamma}
\langle r(t_a)r(t_b) \rangle=&
  \int_{-\infty}^{\infty} \frac{\d f}{(2\pi f)^2}\,(1-e^{-2\pi i f t_a}) (1-e^{2\pi i f t_b})   \sum_{L_1 M_1 L_2 M_2}
       J_{L_1}(x_a) J^*_{L_2}(x_b)\, \times \nonumber \\
       &
       Y_{L_1 M_1}(\hat{e}_a) Y^*_{L_2 M_2}(\hat{e}_b)\,
       \int_{S^2} \d \hat{k} \; R(f,\hat{k},\hat{e}_a,\hat{e}_b)\, Y_{L_1 M_1}^*(\hat{k}) Y_{L_2 M_2}(\hat{k})\,,
\end{align}
\end{widetext}
where $x_a=2\pi f D_a/c$, $x_b=2\pi f D_b/c$, and
\begin{align}
& J_{L_1}(x_a)=  i^{L_1} 2\pi
   \int_0^{x_a}\d x\,e^{ix} j_{L_1}(x),  \\
& J_{L_2}(x_b)=  i^{L_2} 2\pi
   \int_0^{x_b}\d x\,e^{ix} j_{L_2}(x),  \\
&R(f,\hat{k},\hat{e}_a,\hat{e}_b)=
    \sum_{AA'}
    P_{AA'}(f,\hat{k}) \,d_a^{ij} d_b^{kl}\,
    \mathbf{e}^{A}_{ij}(\hat{k}) \mathbf{e}^{*A'}_{kl}(\hat{k})\,.
\end{align}

In terms of the Stokes parameters in Eq.~(\ref{eq:spIQUV}) and the definitions,
\begin{align}
    \mathbb{E}^{I}_{ijkl} (\hat{k}) &=
    \mathbf{e}^{R}_{ij}(\hat{k}) \mathbf{e}^{*R}_{kl}(\hat{k}) +
    \mathbf{e}^{L}_{ij}(\hat{k}) \mathbf{e}^{*L}_{kl}(\hat{k}) \,,
    \nonumber \\
    \mathbb{E}^{V}_{ijkl} (\hat{k}) &=
    \mathbf{e}^{R}_{ij}(\hat{k}) \mathbf{e}^{*R}_{kl}(\hat{k}) -
    \mathbf{e}^{L}_{ij}(\hat{k}) \mathbf{e}^{*L}_{kl}(\hat{k}) \,,
    \nonumber \\
    \mathbb{E}^{Q+iU}_{ijkl} (\hat{k}) &=
    \mathbf{e}^{L}_{ij}(\hat{k}) \mathbf{e}^{*R}_{kl}(\hat{k}) \,,
    \nonumber \\
    \mathbb{E}^{Q-iU}_{ijkl} (\hat{k}) &=
    \mathbf{e}^{R}_{ij}(\hat{k}) \mathbf{e}^{*L}_{kl}(\hat{k}) \,,
\end{align}

denoting $\mathbb{D}$ as the direct product of two detector tensors
\begin{equation}
{}^{ijkl}\mathbb{D}_{ab}=d_a^{ij} d_b^{kl} \,,
\end{equation}
the antenna pattern for signal reception can be rewritten as
\begin{equation}
 R(f,\hat{k},\hat{e}_a,\hat{e}_b)=
\sum_{X=\{I,V,Q\pm iU\}}
 X(f,\hat{k}) \;
{}^{ijkl}\mathbb{D}_{ab} \;
 \mathbb{E}^{I}_{ijkl} (\hat{k}) \,.
\end{equation}
We further expand the Stokes parameters and the polarization basis tensors in terms of ordinary and spin-weighted spherical harmonics as
\begin{align}
    I(f,\hat{k}) &= \sum_{\ell m}I_{\ell m}(f) \; Y_{\ell m}(\hat{k}) \,,\nonumber \\
    V(f,\hat{k}) &= \sum_{\ell m}V_{\ell m}(f) \; Y_{\ell m}(\hat{k}) \,,\nonumber \\
     (Q+iU)(f,\hat{k}) &= \sum_{\ell m}(Q+iU)_{\ell m}(f) \; _{+4}Y_{\ell m}(\hat{k}) \,,\nonumber \\
    (Q-iU)(f,\hat{k}) &= \sum_{\ell m}(Q-iU)_{\ell m}(f) \; _{-4}Y_{\ell m}(\hat{k}) \,,
\end{align}
and
\begin{align}
    \mathbb{E}^{I}_{ijkl} (\hat{k})
    &=\sum_{\ell_e m_e} {}_{ijkl}\mathbb{E}^{I}_{\ell_e m_e} Y_{\ell_e m_e}(\hat{k}) \,, \nonumber
   \\
    \mathbb{E}^{V}_{ijkl} (\hat{k})
    &=\sum_{\ell_e m_e} {}_{ijkl}\mathbb{E}^{V}_{\ell_e m_e} Y_{\ell_e m_e}(\hat{k}) \,, \nonumber
   \\
    \mathbb{E}^{Q+iU}_{ijkl} (\hat{k})
    &=\sum_{\ell_e m_e} {}_{ijkl}\mathbb{E}^{Q+iU}_{\ell_e m_e} {}_{-4}Y_{\ell_e m_e}(\hat{k}) \,, \nonumber
   \\
    \mathbb{E}^{Q-iU}_{ijkl} (\hat{k})
    &=\sum_{\ell_e m_e} {}_{ijkl}\mathbb{E}^{Q-iU}_{\ell_e m_e} {}_{+4}Y_{\ell_e m_e}(\hat{k}) \,,
\end{align}
where the specific combinations, $Q\pm iU$, make them become spin $\pm4$ objects so that we can expand them nicely by the corresponding spin-weighted spherical harmonics. A brief introduction to the spin-weighted spherical harmonics is found in Appendix~\ref{sec:spinweight}. In Appendix~\ref{sec:polartensor}, we have given the multipole moments of the polarization tensors $\mathbb{E}_{ijkl}$ for each Stokes parameter.

Hence we can express the time-residual correlation in the following form,
\begin{widetext}
\begin{equation}
\langle r(t_a)r(t_b) \rangle=
  \int_{-\infty}^{\infty} \frac{\d f}{(2\pi f)^2}\,(1-e^{-2\pi i f t_a}) (1-e^{2\pi i f t_b})
   \sum_{X=\{I,V,Q\pm iU\}}
   \sum_{\ell m}  X_{\ell m}(f) \gamma_{\ell m}^{I}(x_a,x_b,\hat{e}_a,\hat{e}_b)\,,
\end{equation}
where the overlap reduction functions (ORFs) are given by
\begin{align}
    \label{gammaIV_lm}
    \gamma_{\ell m}^{I,V}(x_a,x_b,\hat{e}_a,\hat{e}_b)
  &= \sum_{L_1 M_1 L_2 M_2}
       J_{L_1} J^*_{L_2} Y_{L_1 M_1}(\hat{e}_a) Y^*_{L_2 M_2}(\hat{e}_b)
     \sum_{\ell_e m_e} \mathbb{D} \cdot \mathbb{E}^{I,V}_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && \ell \\
        M_1 && M_2 && 0\; m_e  && 0\;  m
    \end{matrix}
    \right\rangle
    \,, \\
    \label{gammaQU_lm}
 \gamma_{\ell m}^{Q\pm iU}(x_a,x_b,\hat{e}_a,\hat{e}_b)
  &= \sum_{L_1 M_1 L_2 M_2}
       J_{L_1} J^*_{L_2} Y_{L_1 M_1}(\hat{e}_a) Y^*_{L_2 M_2}(\hat{e}_b)
     \sum_{\ell_e m_e} \mathbb{D} \cdot \mathbb{E}^{Q\pm iU}_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && \ell \\
        M_1 && M_2 && \mp4\; m_e  && \pm4\;  m
    \end{matrix}
     \right\rangle
    \,.
\end{align}
\end{widetext}
In Eqs.~(\ref{gammaIV_lm}) and (\ref{gammaQU_lm}), we have replaced the antenna pattern functions by
\begin{equation}
\mathbb{D} \cdot \mathbb{E} \equiv  {}^{ijkl}\mathbb{D}_{ab} \mathbb{E}_{ijkl}\,,
\end{equation}
and introduced a shorthand notation for the integral of a product of four spherical harmonics:
\begin{widetext}
\begin{equation}
    \left\langle
    \begin{matrix}
        L_1 && L_2 && l_1  &&  l_2 \\
        M_1 && M_2 && s_1\;m_1  && s_2\;m_2
    \end{matrix}
    \right\rangle
    \equiv
    \int \d \hat{k} \; Y_{L_1 M_1}^*(\hat{k}) Y_{L_2 M_2}(\hat{k}) \;
     {}_{s_1}\!Y_{l_1 m_1}(\hat{k}) \; {}_{s_2}\!Y_{l_2 m_2}(\hat{k}) \,.
\end{equation}
\end{widetext}
In Appendix~\ref{sec:Wigner3j}, we show that this integral can be expressed as a product of Wigner-3j symbols:
\begin{widetext}
\begin{align}
 \left\langle
    \begin{matrix}
        L_1 && L_2 && \ell_e  &&  \ell \\
        M_1 && M_2 && s\;m_e  && -s\;m
    \end{matrix}
    \right\rangle
=& \sum_{LM}\;
 (-1)^{M_1}\sqrt{\frac{(2L_1+1)(2L_2+1)(2L+1)}{4\pi}}
    \begin{pmatrix}
        L_1 && L_2  &&  L \\
        0 && 0  &&  0
    \end{pmatrix}
    \begin{pmatrix}
        L_1 && L_2  &&  L \\
      -M_1 && M_2  &&  M
    \end{pmatrix}\; \times\nonumber \\
 &   (-1)^M   \sqrt{\frac{(2L+1)(2\ell_e+1)(2\ell+1)}{4\pi}}
    \begin{pmatrix}
        L && \ell_e  &&  \ell \\
        0 && -s  &&  s
    \end{pmatrix}
    \begin{pmatrix}
        L && \ell_e  &&  \ell \\
       -M && m_e  &&  m
    \end{pmatrix} \,.
    \label{4product}
\end{align}
\end{widetext}
Eqs.~(\ref{gammaIV_lm}) and (\ref{gammaQU_lm}) are the most general ORFs for a pair of Galactic pulsars
$a$ and $b$, respectively at distances $D_a$ and $D_b$ from the Earth.

\section{Overlap Reduction Functions in the Computational Frame}

Let us choose the so-called computational frame such that pulsar $a$ is located at the north pole of the Earth
while pulsar $b$ is stayed on the $\phi=0$ meridian. Their coordinates are then given by
\begin{equation}
\hat{e}_a=(0,0,1),\quad \hat{e}_b=(\sin\zeta,0,\cos\zeta)\,
\end{equation}
where $\zeta$ is their separation angle. In this coordinate system, the multipole moments of the antenna pattern functions denoted by $\mathbb{D}_0 \cdot \mathbb{E}$ are listed in Appendix~\ref{sec:DE}. Furthermore, we have
\begin{equation}
Y_{L_1 M_1}(\hat{e}_a)= \sqrt{\frac{2L_1+1}{4\pi}} \delta_{M_1 0},\;\;
Y^*_{L_2 M_2}(\hat{e}_b)=Y^*_{L_2 M_2}(\zeta,0).
\end{equation}
Thus, the integral of the product~(\ref{4product}) simplifies to
\begin{widetext}
\begin{align}
 \left\langle
    \begin{matrix}
        L_1 && L_2 && \ell_e  &&  \ell \\
        M_1 && M_2 && s\;m_e  && -s\;m
    \end{matrix}
    \right\rangle
 =&  \left\langle
    \begin{matrix}
        L_1 && L_2 && \ell_e  &&  \ell \\
        0 && M_2 && s\;m_e  && -s\;m
    \end{matrix}
    \right\rangle
    \nonumber \\
=& \sum_{L}\;
 (-1)^{M_2}\sqrt{\frac{(2L_1+1)(2L_2+1)(2L+1)}{4\pi}}
    \begin{pmatrix}
        L_1 && L_2  &&  L \\
        0 && 0  &&  0
    \end{pmatrix}
    \begin{pmatrix}
        L_1 && L_2  &&  L \\
        0 && M_2  &&  -M_2
    \end{pmatrix}\; \times\nonumber \\
 &    \sqrt{\frac{(2L+1)(2\ell_e+1)(2\ell+1)}{4\pi}}
    \begin{pmatrix}
        L && \ell_e  &&  \ell \\
        0 && -s  &&  s
    \end{pmatrix}
    \begin{pmatrix}
        L && \ell_e  &&  \ell \\
       M_2 && m_e  &&  m
    \end{pmatrix} \,.
\end{align}
\end{widetext}

Hence, Eqs.~(\ref{gammaIV_lm}) and (\ref{gammaQU_lm}) become
\begin{widetext}
\begin{align}
    \label{0gammaIV_lm}
    \gamma_{\ell m}^{I,V}(\zeta)
  &= \sum_{L_1 L_2 M_2}   J_{L_1} J^*_{L_2}
      \sqrt{\frac{2L_1+1}{4\pi}} Y^*_{L_2 M_2}(\zeta,0)
     \sum_{\ell_e m_e} \mathbb{D}_0 \cdot \mathbb{E}^{I,V}_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && \ell \\
        0  && M_2 && 0\; m_e  && 0\;  m
    \end{matrix}
    \right\rangle
    \,, \\
    \label{0gammaQU_lm}
 \gamma_{\ell m}^{Q\pm iU}(\zeta)
  &= \sum_{L_1 L_2 M_2}   J_{L_1} J^*_{L_2}
      \sqrt{\frac{2L_1+1}{4\pi}}  Y^*_{L_2 M_2}(\zeta,0)
     \sum_{\ell_e m_e} \mathbb{D}_0 \cdot \mathbb{E}^{Q\pm iU}_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && \ell \\
        0 && M_2 && \mp4\; m_e  && \pm4\;  m
    \end{matrix}
     \right\rangle
    \,.
\end{align}
\end{widetext}

In Appendix~\ref{sec:DE}, we have found explicit analytic forms for the antenna pattern functions in the computational frame,
$ \mathbb{D}_0 \cdot \mathbb{E}^{I,V}_{\ell_e m_e}$ and $\mathbb{D}_0 \cdot \mathbb{E}^{Q\pm iU}_{\ell_e m_e}$.
Every nonvanishing antenna pattern function is equal to the associated Legendre polynomial $P_2^{m_e} (\cos\zeta)$ times a real constant.
Using Eq.~(\ref{Yconjugate}) and the conjugate relations for the antenna pattern functions in Appendix~\ref{sec:DE}, it is straightforward to show that the four ORFs in the computational frame have the conjugate relations:
\begin{align}
\gamma_{\ell -m}^{I}=&(-1)^{m} \gamma_{\ell m}^{I}\,,
\label{Iconjugate}\\
\gamma_{\ell -m}^{V}=&(-1)^{m+1}\gamma_{\ell m}^{V}\,,
\label{Vconjugate}\\
\gamma_{\ell -m}^{Q\pm iU}=&(-1)^{m} \gamma_{\ell m}^{Q\mp iU}\,.
\label{QUconjugate}
\end{align}
Note that the four ORFs are in general complex functions of $\zeta$. They become real when only the Earth term is considered~\cite{mingar13,kato16}. Our method includes the effects of the pulsar term that give rise to nonzero imaginary parts.

\subsection{Unpolarized Isotropic Case}

In the case of isotropic and unpolarized SGWB, the only relevant ORF is the $\gamma_{00}^{I}$, which can be calculated from Eq.~(\ref{0gammaIV_lm}) as
\begin{widetext}
\begin{align}
   \gamma_{00}^I(\zeta)
  =& \sum_{L_1 L_2 M_2}
      J_{L_1} J^*_{L_2} \sqrt{\frac{2L_1+1}{4\pi}} Y^*_{L_2 M_2}(\zeta,0)
     \sum_{\ell_e m_e} \mathbb{D}_0 \cdot \mathbb{E}^I_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && 0 \\
        0  && M_2 && 0\; m_e  && 0\;  0
    \end{matrix}
    \right\rangle
     \nonumber \\
   =& \sum_{L_1 L_2 \ell_e m_e}  J_{L_1} J^*_{L_2}
      (-1)^{m_e} \frac{1}{4\pi}  Y_{L_2 m_e}(\zeta,0)\;
      \mathbb{D}_0 \cdot \mathbb{E}^I_{\ell_e m_e} \times \nonumber \\
     &\sqrt{\frac{(2L_1+1)^2(2L_2+1)(2\ell_e+1)}{4\pi}}
    \begin{pmatrix}
        L_1 && L_2  &&  \ell_e \\
        0 && 0  &&  0
    \end{pmatrix}
      \begin{pmatrix}
        L_1 && L_2  &&  \ell_e \\
       0  && -m_e  &&  m_e
    \end{pmatrix} \,.
\end{align}
\end{widetext}
We have numerically computed $\gamma_{00}^I(\zeta)$ using $x_a=x_b=20\pi$ and summing $L_2$ up to $80$. Note that for each value of $L_2$, the range of $L_1$ is determined by the triangular condition in the Wigner-3j symbol. Figure~\ref{fig:gammaI00} shows $\gamma_{00}^I(\zeta)$ plotted against the separation angle $\zeta$. At large angular separation the real part of $\gamma_{00}^I(\zeta)$ reproduces the Hellings and Downs curve
for the quadrupolar interpulsar correlations~\cite{downs},
which has considered the Earth term only, while the pulsar term contributes to a relatively small imaginary part.
$\gamma_{00}^I(\zeta)$ gains power from the pulsar term and begins to deviate from the Hellings and Downs curve at small angular separation. The autocorrelation $\gamma_{00}^I(\zeta=0^\circ)$ has power two times larger than the Hellings and Downs curve. This small-angle behavior is consistent with the results by considering the corrections from the pulsar term~\cite{mingar14}.

\subsection{Circularly Polarized Isotropic Case}
\label{CPI}

In the case of isotropic and circularly polarized SGWB, the other relevant ORF is the $\gamma_{00}^{V}$, which can be calculated from Eq.~(\ref{0gammaIV_lm}) as
\begin{widetext}
\begin{align}
   \gamma_{00}^V(\zeta)
  =& \sum_{L_1 L_2 M_2}  J_{L_1} J^*_{L_2}
       \sqrt{\frac{2L_1+1}{4\pi}} Y^*_{L_2 M_2}(\zeta,0)
     \sum_{\ell_e m_e} \mathbb{D}_0 \cdot \mathbb{E}^V_{\ell_e m_e}
    \left\langle
    \begin{matrix}
        L_1 && L_2  && \ell_e   && 0 \\
        0  && M_2 && 0\; m_e  && 0\;  0
    \end{matrix}
    \right\rangle
     \nonumber \\
    =& \sum_{L_1 L_2 \ell_e m_e} J_{L_1} J^*_{L_2}
      (-1)^{m_e} \frac{1}{4\pi}  Y_{L_2 m_e}(\zeta,0)\;
      \mathbb{D}_0 \cdot \mathbb{E}^V_{\ell_e m_e} \times \nonumber \\
     & \sqrt{\frac{(2L_1+1)^2(2L_2+1)(2\ell_e+1)}{4\pi}}
    \begin{pmatrix}
        L_1 && L_2  &&  \ell_e \\
        0 && 0  &&  0
    \end{pmatrix}
      \begin{pmatrix}
        L_1 && L_2  &&  \ell_e \\
       0  && -m_e  &&  m_e
    \end{pmatrix} \,.
\end{align}
\end{widetext}
We will show that $\gamma_{00}^V(\zeta)=0$ as follows. There is no contribution to the summation from terms with $m_e=0$ because $\mathbb{D}_0 \cdot \mathbb{E}^V_{\ell_e 0}=0$. Furthermore, since we have
$\mathbb{D}_0 \cdot \mathbb{E}^V_{\ell_e -m_e}=(-1)^{m_e+1} \mathbb{D}_0 \cdot \mathbb{E}^V_{\ell_e m_e}$
and the property,
\begin{equation}
 \begin{pmatrix}
        l_1 && l_2  &&  l_3 \\
       -m_1  && -m_2  && - m_3
  \end{pmatrix}
  = (-1)^{l_1 + l_2 +  l_3}
  \begin{pmatrix}
        l_1 && l_2  &&  l_3 \\
       m_1  && m_2  &&  m_3
  \end{pmatrix} \,,
\end{equation}
which is zero when $m_1=m_2=m_3=0$ unless $l_1 + l_2 +  l_3$ is an even integer, the contribution from terms with positive values of $m_e$ exactly cancels that from terms with negative values of $m_e$,
making the overall summation vanish.

\subsection{Higher Multipole Moments}
\begin{widetext}

\begin{figure}[t]
\centering
\includegraphics[width=0.48\textwidth]{gamma00I_SmallAngle.pdf}
\includegraphics[width=0.48\textwidth]{gamma00I_LargeAngle.pdf}
\caption[Coordinate]{The right panel shows $\gamma_{00}^I(\zeta)$, while the left panel shows the power at small angular separation. We have used $x_a=x_b=20\pi$. The dashed curve is the Hellings and Downs curve. The real part of
$\gamma_{00}^I(\zeta)$ (solid curves) overlaps with the Hellings and Downs curve at large angular separation, while the imaginary part of $\gamma_{00}^I(\zeta)$ (dotted curves) receives a small contribution from the pulsar term.}
\label{fig:gammaI00}
\index{figures}
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.48\textwidth]{gamma1mI.pdf}
\includegraphics[width=0.48\textwidth]{gamma2mI.pdf}
\caption[Coordinate]{$\gamma_{\ell m}^I(\zeta)$ for $\ell=1,2$ and $m\ge 0$, with $x_a=x_b=20\pi$. The multipoles with $m<0$ are given by the relation~(\ref{Iconjugate}), $\gamma_{\ell -m}^{I}=(-1)^{m} \gamma_{\ell m}^{I}$.}
\label{fig:gammaIlm}
\index{figures}
\end{figure}
\end{widetext}
We have numerically computed $\gamma_{\ell m}^I(\zeta)$ using $x_a=x_b=20\pi$ and summing $L_2$ up to $8$ for
$\ell=1,2$, as shown in Fig.~\ref{fig:gammaIlm}. Note that we have plotted the ORFs for $\zeta>5^\circ$. To obtain the small-angle resolution that shows the effects from the pulsar term, we would need to increase the maximum value of $L_2$. Similarly, Fig.~\ref{fig:gammaVlm} shows $\gamma_{\ell m}^V(\zeta)$ for $\ell=1,2$. Note that $\gamma_{\ell 0}^V(\zeta)=0$, which can be shown following the same reasoning as in Sec.~\ref{CPI}. These ORF harmonics reproduce those found in Refs.~\cite{mingar13,kato16} using the Shapiro time delay with the Earth term only. Figure~\ref{fig:gammaQUlm} shows our new results for $\gamma_{\ell m}^{Q\pm iU}(\zeta)$ for $\ell=4$. In these figures, we have omitted the imaginary parts that are small compared to the real parts.
\begin{widetext}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.48\textwidth]{gamma1mV.pdf}
\includegraphics[width=0.48\textwidth]{gamma2mV.pdf}
\caption[Coordinate]{$\gamma_{\ell m}^V(\zeta)$ for $\ell=1,2$ and $m\ge 0$, with $x_a=x_b=20\pi$. The multipoles with $m<0$ are given by the relation~(\ref{Vconjugate}), $\gamma_{\ell -m}^{V}=(-1)^{m+1} \gamma_{\ell m}^{V}$.}
\label{fig:gammaVlm}
\index{figures}
\end{figure}

\begin{figure}[ht!]
\centering
\includegraphics[width=0.48\textwidth]{gammalmQpiU.pdf}
\includegraphics[width=0.48\textwidth]{gammalmQmiU.pdf}
\caption[Coordinate]{$\gamma_{\ell m}^{Q\pm iU}(\zeta)$ for $\ell=4$ and $m\ge 0$, with $x_a=x_b=20\pi$.
The multipoles with $m<0$ are given by the relation~(\ref{QUconjugate}),
$\gamma_{\ell -m}^{Q\pm iU}=(-1)^{m} \gamma_{\ell m}^{Q\mp iU}$.}
\label{fig:gammaQUlm}
\index{figures}
\end{figure}
\end{widetext}

\section{Overlap Reduction Functions in the Celestial Coordinates}

We can readily get the ORFs in the celestial coordinates from those in the computational frame through a three-dimensional rotation~\cite{chu21}. Let pulsar $a$ be located at the polar angles $\hat{e}_a(\theta,\phi)$, and the angle between its longitude and the great circle connecting the pulsar pair $a$-$b$ be $\alpha$. Then, Eqs.~(\ref{gammaIV_lm}) and (\ref{gammaQU_lm}) are connected to Eqs.~(\ref{0gammaIV_lm}) and (\ref{0gammaQU_lm}) by
\begin{align}
    \gamma_{\ell m}^{I,V}(\hat{e}_a,\hat{e}_b)
  &= \sum_{m'} D^{\ell}_{m' m}(-\alpha,-\theta,-\phi) \gamma_{\ell m'}^{I,V}(\zeta)   \,, \\
 \gamma_{\ell m}^{Q\pm iU}(\hat{e}_a,\hat{e}_b)
  &=  \sum_{m'} D^{\ell}_{m' m}(-\alpha,-\theta,-\phi) \gamma_{\ell m'}^{Q\pm iU}(\zeta)   \,,
\end{align}
where the three-dimensional rotation is represented by the Wigner-$D$ matrix, which is closely related to the spin-weighted spherical harmonics:
\begin{align}
    \label{eq:DtoSWSH}
    D^{\ell}_{s m}(\alpha,\beta,\gamma)
    = \sqrt{\frac{4\pi}{2\ell+1}} {}_{-s}Y_{\ell m}(-\beta,-\gamma) e^{-i s \alpha} \,.
\end{align}

\section{Conclusion}

We have developed a numerical scheme to compute the ORFs for the Stokes parameters of an anisotropic stochastic gravitational wave background in pulsar-timing-array observation. We have used the Sachs-Wolfe line-of-sight integral to construct the time-residual correlation between a pair of Galactic millisecond pulsars on the sky. The integral incorporates simultaneously the effects of the pulsar term that contribute power to the ORFs at small angular separation of the pulsar pair. The ORFs for the linear polarization are given for the first time. Based on a spherical harmonic analysis, our method allows us to compute the ORF multipoles order by order as well as improving their resolution at small angular separation in a controllable manner. The method can be readily applied to compute the ORFs for a Galactic pulsar pair at different distances from the Earth. Furthermore, the Sachs-Wolfe line-of-sight integral can be used to consider extragalactic pulsars at high redshifts.

\begin{acknowledgments}
This work was supported in part by the Ministry of Science and Technology (MOST) of Taiwan, Republic of China, under
Grants No. MOST 109-2112-M-032-006 (G.C.L.) and No. MOST 109-2112-M-001-003 (K.W.N.).
\end{acknowledgments}