\section*{Abstract}
\subsection*{Background}
Clinical research and drug development in orphan diseases is challenging, since large-scale randomized studies are difficult to conduct. Formally synthesizing the evidence is therefore of great value, yet this is rarely done in the drug approval process. Phase III designs that make better use of phase II data can facilitate drug development in orphan diseases.

\subsection*{Methods}
A Bayesian meta-analytic approach is used to inform the phase III study with phase II data. It is particularly attractive, since uncertainty of between-trial heterogeneity can be dealt with probabilistically, which is critical if the number of studies is small. Furthermore, it allows quantifying and discounting the phase II data through the predictive distribution relevant for phase III. A phase III design is proposed which uses the phase II data and considers approval based on a phase III interim analysis. The design is illustrated with a non-inferiority case study from an FDA approval in herpetic keratitis (an orphan disease). Design operating characteristics are compared to those of a traditional design, which ignores the phase II data.

\subsection*{Results}
An analysis of the phase II data reveals good but insufficient evidence for non-inferiority, highlighting the need for a phase III study. For the phase III study supported by phase II data, the interim analysis is based on half of the patients. For this design, the meta-analytic interim results are conclusive and would justify approval. In contrast, based on the phase III data only, interim results are inconclusive and would require further evidence.

\subsection*{Conclusions}
To accelerate drug development for orphan diseases, innovative study designs and appropriate methodology are needed. Taking advantage of randomized phase II data when analyzing phase III studies looks promising because the evidence from phase II supports informed decision making. The implementation of the Bayesian  design is straightforward with public software such as \texttt{R}.\\[2\baselineskip]
\textbf{Keywords:} drug development in rare diseases; phase III studies; Bayesian statistics; meta-analysis
\clearpage

\section*{Introduction}
\label{sec:introduction}
\input{introduction}

\section*{Methods}
\label{sec:methods}
\input{methods}

\section*{Case Study}
\label{sec:application}
\input{application}

\section*{Discussion}
\label{sec:discussion}
\input{discussion}

\clearpage
\section*{Declaration of Conflicting Interest}
Dr. Wandel and Dr. Neuenschwander are employed by Novartis Pharma AG, Basel, Switzerland.

\clearpage
\input{bibliography}

\clearpage
\input{tables}

\clearpage
\input{figures}

\clearpage
\input{appendix}
\end{document}

\section*{Appendix}
\begin{verbatim}
library("bayesmeta")
library("metafor")

################################################
#  Part I - ratio of RRs for tau ~ HN(...)
################################################
set.seed(314)
N <- 1000000
tau1 <- abs(rnorm(N, 0, 0.5))
tau2 <- abs(rnorm(N, 0, 1.0))
theta1 <- rnorm(N, 0, tau1)
theta2 <- rnorm(N, 0, tau2)
exp(quantile(theta1, prob = 0.975))
exp(quantile(theta2, prob = 0.975))

################################################
#  Part II - application
################################################

# --------------------------------------------
# read-in data
# transform into 2x2 cell entries
# derive mean and se from normal approximation
# --------------------------------------------

all.data <- data.frame(study = c("4", "5", "6", "7IA", "7FA"),
                       rt    = c(19, 15, 31, 35, 74),
                       nt    = c(23, 18, 36, 40, 84),
                       rc    = c(16, 12, 27, 36, 73),
                       nc    = c(22, 17, 38, 40, 80))

all.data$ai <- all.data$rt
all.data$bi <- all.data$nt - all.data$rt
all.data$ci <- all.data$rc
all.data$di <- all.data$nc - all.data$rc

nmappr <- escalc(ai=ai, bi=bi, ci=ci, di=di,
                 data=all.data,	measure="RR")

# --------------------------------------------------
# analyses:
#  - end-of-phase-II
#  - interim for phase III
#  - final for phase III
# summaries:
#  - population mean,
#  - predicted effect for end-of-phase-II
#  - phase-III effect for interim and final analyses
# --------------------------------------------------

eop2  <- bayesmeta(y = nmappr$yi[1:3], sigma = sqrt(nmappr$vi[1:3]),
                   tau.prior = function(t){dhalfnormal(t,scale=0.5)})

ph3IA <- bayesmeta(y = nmappr$yi[1:4], sigma = sqrt(nmappr$vi[1:4]),
                   tau.prior = function(t){dhalfnormal(t,scale=0.5)})

ph3FA <- bayesmeta(y = nmappr$yi[c(1:3,5)], sigma = sqrt(nmappr$vi[c(1:3,5)]),
                   tau.prior = function(t){dhalfnormal(t,scale=0.5)})

round(exp(t(eop2$summary)[,c("median", "9
round(exp(ph3IA$theta[c("median", "9
round(exp(ph3FA$theta[c("median", "9
\end{verbatim}
\end{appendix}

\subsection*{Background}
Herpetic keratitis is an inflammatory condition of the eye caused by an outbreak of the herpes simplex virus (HSV)\cite{Kay2006,White2014}. It can have serious consequences and remains the leading cause of corneal blindness in the industrialized world \cite{Dawson1976, Suresh1999}. With as few as 1.5 million people affected world-wide~\cite{Farooq2012}, it has been classified as an orphan indication by the FDA~\cite{CorrespFDA} and the European Medical Agency~\cite{EMAOrphan2008}. \\[2\baselineskip]

In 2009, the FDA approved Zirgan for the treatment of herpetic keratitis (dendritic ulcers)~\cite{ApprovalFDA2009}. To discuss all details of the approval is beyond the scope of this application (see the publicly available documents~\cite{ApprovalDocsFDA}). However, a few points are noteworthy. Most importantly, from the files~\cite{CorrespFDA,ApprovalDocsFDA} it appears that approval was based on a retrospective analysis of the four relevant studies, three phase II and one phase III study. Retrospective means that the sponsor submitted the results of the studies after they were conducted, rather than seeking the agency's advice beforehand. Subsequently, this led to discrepancies between the sponsor's and FDA's primary analyses, including changes of the population, of the endpoint and from superiority to non-inferiority.\\[2\baselineskip]

The reasons behind this rather unusual approach to approval are not entirely clear. One explanation may be that the original manufacturer (Th\'{e}a of France) did not intend to bring Zirgan to the US market on its own; rather, it sold the license for the US market to Sirion Therapeutics in 2007 which then initiated the submission. This and the fact that the clinical studies were already conducted in the 1990s may explain why no early discussions with the FDA took place.\\[2\baselineskip]

Our goal here is not to reconstruct the approval history in detail. Rather, we will use the example to discuss an alternative, more efficient statistical approach towards approval, based on the following design specifications in the non-inferiority setting: cure rate at day 14 as endpoint, dendritic and geographic ulcers as population, and an absolute non-inferiority margin of 12 percentage points. Furthermore, we will use the risk ratio (RR) to quantify the treatment effect.\\[2\baselineskip]

In the following, we present the evidence available at the hypothetical end-of-phase II meeting, a potential phase III trial and approval strategy, and the results of the actual phase III trial.

\subsection*{Hypothetical end-of-phase II meeting}
Three randomized phase II studies~\cite{ZirganClinPaper} were conducted between April 1990 and October 1992 (Table~\ref{tbl:data}). The studies were similar, with the only minor difference being the treatment regimen in study 6. For simplicity, we assume that this difference is not relevant for the clinical outcome.\\[2\baselineskip]

We now turn the clock back and assume we are in the situation of an end-of-phase II meeting. We assume that the sponsor would agree to a non-inferiority analysis of Zirgan versus Acyclovir (the standard of care) with the primary endpoint being cure rate at day 14. Actually, setting a non-inferiority margin proved to be difficult. For cure rate at day 14, the FDA determined two effect sizes M1~\cite{FDAGuidanceNI}: 14\% and 18\%~\cite{ApprovalDocsFDA}. The latter implies an absolute non-inferiority margin of 12 percentage points when retaining one third of the effect. We assume here that this margin had been agreed to.\\[2\baselineskip]

At this stage, it is interesting to perform a non-inferiority analysis (Zirgan versus Acyclovir) of the phase II data. If the evidence were overwhelming, it would be fair to ask whether a phase III study were required, or if approval could be granted based on the phase II data only.\\[2\baselineskip]

Our interest is the phase III treatment effect. However, since no phase III data are available yet, the phase III treatment effect corresponds to the predicted treatment effect $\theta_{\star}$ from the phase II studies (see Section 2). The underlying statistical model is the NNHM~\eqref{EX},~\eqref{EX:normal}, with study-specific estimates of the log-risk-ratios $Y_j = \log(\text{RR}_j)$ and standard errors

\begin{equation}
s_j = \sqrt{1/r_C - 1/n_C + 1/r_T - 1/n_T}
\end{equation}

 where $n$ and $r$ denote the number of patients and responders. This requires a transformation to the risk difference scale and a sensible prior distribution for the between-trial heterogeneity parameter $\tau$ (the prior for $\mu$ will be non-informative).\\[2\baselineskip]

The first point is straightforward. For a response rate $p_C$ in the control group and a pre-defined non-inferiority margin $m=p_C - p_T$, the transformation is given by the definition of the risk ratio; for $p_C = 0.9$ (the assumed cure rate for Acyclovir based on historical data) and margin $m = 0.12$, non-inferiority holds if $\text{RR}_{T:C} \geq 0.867$.\\[2\baselineskip]

For the $\tau$ prior we use $\tau \sim HN(0.5)$, which has median 0.34 and 95\% interval (0.016;1.12). This prior is centered at moderate to substantial heterogeneity and covers small to large heterogeneity (see Table~\ref{tab:heterogeneity}). Notably, one may perform a sensitivity analysis using a prior which favors larger between-trial heterogeneity, e.g. $\tau \sim HN(1)$.\\[2\baselineskip]

The meta-analysis of the phase II data is shown in Figure~\ref{fig:eop2_ma}, where the data, study-specific (stratified) risk ratios $\text{RR}_j$, the  population mean $\mu$ and the predicted effect $\theta_{\star}$ are shown. The reference line is drawn at the non-inferiority margin ($\text{RR}_{T:C}=0.867$ for $p_C = 0.9$). The posterior for $\tau$ indicates small between-trial heterogeneity, with median 0.12 (95\% interval 0.00 to 0.51). \\[2\baselineskip]

The meta-analysis provides evidence for non-inferiority. If $\mu$ were the parameter of interest, an almost conclusive statement would follow: the lower bound of the 95\% interval is just below the non-inferiority margin. In fact, $P(\mu \geq \log(0.867)) = 97.1\%$, very close to $97.5\%$. However, the parameter $\theta_{\star}$ in the phase III trial is of interest. For this parameter, the evidence for non-inferiority is weaker, but still substantial: $P(\theta_{\star} \geq \log(0.867)) = 92.0\%$.

\subsection*{Phase III study and proposed strategy for approval}
Designing a phase III study that allows to assess non-inferiority in combination with the available evidence is desirable. Not only will this allow to run a smaller study, it will also provide a treatment effect estimate based on all relevant evidence. However, regulators may have good reasons to argue that a smaller study may provide insufficient information for approval, especially to assess the safety and risk/benefit ratio.\\[2\baselineskip]

We now discuss the design of  a phase III study (study 7) which uses phase II data and allows for seeking approval based on an interim analysis. Depending on negotiations with regulators, a post-approval commitment to run the study to its end (even if approval is granted at interim) may be required. However, such negotiations will always be case-specific, highlighting the importance of early discussions with regulators. Nevertheless, the option to seek approval based on a positive interim analysis seems attractive for this case study. Since the endpoint is evaluated at day 14, there will be a small time window between the last patient enrolled for the interim analysis and the actual data read-out and analysis. With an anticipated recruitment period of two years, such a strategy could result in a markedly earlier approval. \\[2\baselineskip]

When seeking approval based on interim results, the information fraction for the interim analysis becomes a key design aspect. We will assume that the interim analysis is conducted after 50\% of the patients have been evaluated. For the sample size, in order to align with the actual study as originally conducted, we will assume $n_C = n_T = 80$. This results in interim sample sizes $\tilde{n}_{C} = \tilde{n}_{T} = 40$.\\[2\baselineskip]

It is also important to understand how much phase II information is borrowed (which depends on the between-trial heterogeneity) when inferring the phase III effect. Using the variance ratio approach (Section 2), the ESS is 14.

\subsection*{Operating characteristics}
We evaluate the operating characteristics (type I error rate and power) of the design and compare them to a phase III design ignoring the phase II data. The operating characteristics are presented in Table~\ref{tbl:oc} based on 10'000 simulations conducted in \texttt{R}~\cite{Rsoftware} with the package \texttt{bayesmeta}~\cite{bayesmeta}. For different response rates $p_C$ and treatment differences $\delta$, two probabilities are shown: the probability to be successful at the final analysis (regardless of the outcome of the interim analysis), and the probability to be successful both at the interim and the final analysis.\\[2\baselineskip]

The gain in power for the proposed design can be substantial. For example, for $p_C = 0.7$ (the Acyclorivr cure rate observed in phase II) and $\delta = 0.06$, the power is 87\% versus 66\%. The power gain is even larger at interim (70\% versus 35\%). When $p_C = 0.9$ (the observed cure rate for Acyclovir based on historical data) and $\delta = 0$, the power is 87\% versus 79\%, and 68\% versus 48\% at interim. The larger increase in power at interim is remarkable and due to the highly consistent phase II results, which suggested superiority of Zirgan.\\[2\baselineskip]

The gain in power, however, comes at the price of an increased type I error rate. Strict type I error rate control cannot be guaranteed~\cite{Viele2014}. For example, for $p_C = 0.7$ and $p_C = 0.9$, the type I error rates are 6\% versus 1\% and 8\% versus 3\%. This increase is not dramatic, yet it cannot be ignored and needs to be discussed with regulators during the design phase. If it is of concern, robust prior distributions could be considered~\cite{Schmidli2014}.

\subsection*{Actual phase III data and analysis}
The actual data observed in the phase III study are only available for the final analysis. In order to reconstruct an interim analysis using half of the patients, we use an interim sample size of 40 per arm. Furthermore, we choose the number of responders such that observed response rate at interim is close to the observed response rate at the final analysis (see Figure~\ref{fig:final_analysis}).\\[2\baselineskip]

The results are presented in Figure~\ref{fig:final_analysis}. The interim analysis based on all data (meta-analysis) allows to declare non-inferiority. Note that non-inferiority is claimed based on the parameter corresponding to study 7 ($\theta_{\star}$) incorporating the evidence from studies 4, 5 and 6. On the other hand, the evidence from the phase III study alone is insufficient to declare non-inferiority at interim because the 95\% interval includes the non-inferiority margin.\\[2\baselineskip]

As mentioned before, the idea would be to gain approval with the interim phase III data supported by phase II via the meta-analysis, assuming other data (such as safety) is also favorable. Yet, depending on negotiations with regulators, the study may still run to its end, allowing a more robust evaluation of the effect at the final analysis. The results for the final analysis are also shown in Figure~\ref{fig:final_analysis}. For the meta-analysis, the interval for the risk ratio becomes narrower and still excludes the non-inferiority margin, thus confirming the interim result. The analysis using the phase III study leads to a lower bound of the interval ($0.870$) which is just above the non-inferiority threshold $0.867$, also allowing to conclude non-inferiority.\\[2\baselineskip]

Finally, results for $\tau$ indicate small between-trial heterogeneity at the interim and the final analysis. The posterior median (95\% interval) is 0.12 (0.00 to 0.41) for the interim and  0.13 (0.00 to 0.43) for the final analysis. This supports the consistency of the results across all studies.\\[2\baselineskip]

\subsection*{Hierarchical models}
Hierarchical models (HM) are widely used when data are available from
more than one trial.  The models have two components: a data model and
a parameter model. Data $Y_j$ from trial $j=1,\ldots,J$ follow a
distribution $F$ parameterized by trial-specific parameters
$\theta_j$
\begin{equation}
\label{data:model}
Y_j \vert \theta_j \sim F(\theta_j)
\end{equation}
and trial parameters $\theta_j$ follow a distribution $G$
\begin{equation}
\label{parameter:model}
\theta_j \vert \eta \sim G(\eta)
\end{equation}
Inference for trial parameters can be done in a classical or Bayesian
way. The simplest hierarchical model assumes (approximately) normal
data. Often, the $Y_j$ are parameter estimates rather
than individual data. For this case, the normal-normal hierarchical
model (NNHM) is widely used:
\begin{equation}
\label{EX}
Y_j \vert \theta_j \sim N(\theta_j,s_j^2)
\end{equation}
and
\begin{equation}
\label{EX:normal}
\theta_j \vert \mu,\tau \sim N(\mu,\tau^2)
\end{equation}
For fixed standard errors $s_j$ and known (assumed) $\tau$, classical
and Bayesian conclusions for $\mu$ and trial parameters $\theta_j$ are
the same if a non-informative (improper) prior for $\mu$ is used. For
precision (inverse-variance) weights $w_j$, total precision $w_+$, and
shrinkage parameters $B_j$
\begin{equation}
w_j = \frac{1}{s_j^2+\tau^2}, \qquad w_+ = \sum_{j=1}^Jw_j,  \qquad B_j = \frac{s_j^2}{s_j^2+\tau^2}
\end{equation}
the posterior distribution of $\mu$ based on $Y=(Y_1,\ldots,Y_J)$ is
\begin{equation}
\label{NNHM:mu}
\mu \vert Y, \tau \sim N(\sum w_jY_j/w_+,1/w_+)
\end{equation}
The posterior distributions of the trial parameters $\theta_j$ are
\begin{equation}
\label{NNHM:theta}
\theta_j \vert Y,\tau \sim N(B_j\hat{\mu}+(1-B_j)Y_j,
B_j(\tau^2+B_j/w_+))
\end{equation}
where $\hat{\mu}$ is the posterior mean in (\ref{NNHM:mu}). Classical
analogues to the posterior means and standard deviations are
maximum-likelihood estimates and their standard errors. The special
cases of complete pooling and stratification arise for $\tau=0$ and
$\tau=\infty$, respectively.\\[2\baselineskip]

Intermediate values of $\tau$ lead to
different degrees of information sharing across trials, with
the desirable properties one expects from
an approach aiming to improve inference by
borrowing information from similar trials:
\begin{itemize}
\item The hierarchical model shrinks the trial estimates towards the
  estimate of $\mu$, which acts as a
  safeguard against over-interpreting extreme (good or bad) trial results.
  Shrinkage depends on trial size and between-trial heterogeneity. For
  large trials (small
  $s_j$), shrinkage is small, and notable
  shrinkage is only possible if $\tau$ is of small to moderate size.
\item
The hierarchical model improves precision. Since
\begin{equation}
\label{HM:var}
B_j(\tau^2+\frac{B_j}{w_+}) =
s_j^2 - s_j^2B_j(1-\frac{w_j}{w_+})
\end{equation}
the variance in (\ref{NNHM:theta}) is always smaller than the variance $s_j^2$ of $Y_j$.
\end{itemize}

\subsection*{Between-trial heterogeneity}
The degree of between-trial heterogeneity (standard
deviation $\tau$ in (\ref{EX:normal})) for
the parameters $\theta_1,\ldots,\theta_J$ depends on the parameter
scale and the outcome standard deviation $\sigma$ for one observation
unit (for example one subject or one event). Table
\ref{tab:heterogeneity} shows four typical heterogeneities and
respective $\tau$ values for $\sigma=2$, which is
often used as the reference standard deviation for
normal approximations of binomial, count, and survival data~\cite{Spiegelhalter2004}.
For the four heterogeneities, Table
\ref{tab:heterogeneity} shows the range of parameter values expressed
as the ratio between the 97.5\% quantile and the median; for example,
$\tau=1$ implies a ratio of 7.1, which is clearly large and will be
rare in practice. \\[2\baselineskip]

For the common case of few trials, the size of between-trial
heterogeneity is usually highly uncertain because $\tau$ cannot be
inferred well from the data.  Therefore, it is important to use prior
distributions covering plausible $\tau$ values. Half-normal,
half-Cauchy, and half-t distributions have been suggested in this
context~\cite{Spiegelhalter2004, Gelman2006, Polson2012}.  For the
log-risk ratios used in the application, we will consider
half-normal distributions~\cite{Spiegelhalter2004} with scale
parameters 0.5 and 1, which have medians (95\%-intervals) equal to
0.34 (0.016,1.12) and 0.67 (0.031,2.24), respectively. Since $\tau = 1$ represents large heterogeneity, both priors
are weakly informative, covering small to large
heterogeneity and leaving small probabilities to unrealistically large
heterogeneities, whereby the latter prior (with median 0.67) is rather
conservative. For these priors, the 97.5\% quantile to median ratio for
risk ratio (RR) trial parameters is 2.98 and 8.89, respectively (see Appendix).

\subsection*{Meta-analytic-predictive (MAP) prior }
When designing a new trial with parameter   $\theta_{\star}$, the predictive distribution based
on previous data   $Y_1,\ldots,Y_J$ constitutes the prior distribution for the new
  trial. This is known as \emph{the meta-analytic-predictive (MAP)}
  prior~\cite{Spiegelhalter2004,Neuenschwander2010,Schmidli2014}
\begin{equation}
\label{map:prior}
\theta_{\star} \vert Y_1,\ldots,Y_J
\end{equation}
For the NNHM with known $\tau$
\begin{equation}
\label{NNHM:thetastar}
\theta_{\star} \vert Y_1,\ldots,Y_J, \tau \sim N(\hat{\mu},\tau^2+1/w_+)
\end{equation}
which follows from (\ref{NNHM:theta}) by adding the new trial (with no
data) to the model, i.e., $s_{\star}=\infty$ and $B_{\star}=1$.

\subsection*{Analysis for new trial }
Eventually, after the new data $Y_{\star}$ have been
observed, inference for $\theta_{\star}$ can be done in two ways:
\begin{itemize}[leftmargin=11.0mm]
\item[\MAP] the \emph{meta-analytic-predictive (MAP)} approach
  formally combines the prior (\ref{map:prior}) with
  $\Ystar$ in a standard Bayesian way.
\item[\MAC] the \emph{meta-analytic-combined (MAC)} approach does not require a prior
distribution for $\thstar$. It simply infers $\thstar$ at the end of the new
trial by a meta-analysis of historical and new data, resulting in
\begin{equation}
\label{mac:prior}
\thstar \vert Y_1,\ldots,Y_J,\Ystar
\end{equation}
\end{itemize}
Importantly, \emph{MAC} and \emph{MAC} give identical
  results~\cite{Schmidli2014}.  The \emph{MAP} approach is technically more involved
  because \emph{MAP} priors (\ref{map:prior}) do not follow standard distributions and are typically heavy-tailed.
  This complicates the Bayesian analysis with $Y_{\star}$ at the end
  of the trial, which can be addressed via mixture approximations~\cite{Schmidli2014}. However, even if a \emph{MAC} analysis will
  usually be the method of choice and easy to perform with
  meta-analytic software, the \emph{MAP} prior plays
  an important role: it quantifies prior information at the design
  stage, which may be required in the trial protocol.

\subsection*{Effective sample sizes}
\label{sec:methods:ess}
In many applications, the appropriate use of prior information will
lead to smaller trials. The amount of information is ideally expressed
as an equivalent approximate prior \emph{effective sample size (ESS)}.
In our setting we are interested in $ESS_{\star}$,
the prior effective sample size of the MAP prior
(\ref{map:prior}). Various approaches to $ESS$ have been
proposed~\cite{Neuenschwander2010, Malec2001, Pennello2008, Morita2008, Hampson2014}; they are similar in the sense that they relate
the $ESS$ to the precision (inverse of variance) of the prior
distribution.\\[2\baselineskip]
Here, we will use an approximate two-variances approach which requires: the
variance $V_{\star}$ of the analysis of interest, for which the
$ESS_{\star}$ is unknown; and, the variance $V_0$ of a simpler
analysis (e.g. a meta-analysis with $\tau=0$) with
known $ESS_0$. Assuming that effective sample sizes are approximately
proportional to precisions, the ESS of interest is
\begin{equation}
ESS_{\star} = ESS_0 \times \frac{V_0}{V_{\star}}
\end{equation}
In our case, $V_{\star}$ will be variance of the $MAP$ prior
(\ref{map:prior}), whereas $V_0$ will be the one from
the analysis assuming no between-trial heterogeneity ($\tau=0$).\\[2\baselineskip]


\section{Introduction}
\label{introduction}

This contribution constitutes but a small part of the two-day splinter session entitled ``Star clusters from space, from the ground, and over time'' which took place at the 19th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun in Uppsala, Sweden in June 2016. Given the nature of this contribution it will be far from comprehensive, both in terms of the range of fundamental stellar parameters covered, but also in terms of the various methods employed to estimate such parameters (for more information on these the reader is referred to other Cool Stars 19 proceedings contributions). A recent and more comprehensive discussion on star clusters which not only covers global properties of both young and old clusters, but also our current understanding of the dynamical evolution of clusters, can be found in the proceedings of the Ecole Evry Schatzman 2015 (EES2015) school ``Stellar clusters: benchmarks of stellar physics and galactic evolution''.

Clusters have long represented benchmarks with regard to the determination of fundamental stellar parameters, in large part due to the underlying assumption that members within such ensembles share several common properties; namely they are coeval, have the same chemical composition and are located at roughly the same distance. In addition to these shared characteristics, it is also observationally advantageous to focus on clusters as they have a significantly higher stellar number density (per unit area on the sky) compared to either young associations/moving groups or field stars and so for a given allocation of telescope time one can thus maximise the number of stars in ones sample.

Studies of clusters have also been instrumental in driving our understanding of the formation and evolution of stars. By studying a given cluster we can infer the mass dependence of astrophysical phenomena at a given epoch and by studying several clusters spanning a range of ages we can track how such phenomena evolve with time, as well as investigate second-order effects such as the local environment.

\section{Global parameters}
\label{cmd}

\begin{figure}[h]
\centering
\includegraphics[width=\columnwidth]{Figs/pleiades_cmd_label.pdf}
 \caption{The $V, V$--$I_{\rm{c}}$ CMD of the Pleiades with several regions marked and the global parameters one can estimate using stars in these regions.}
   \label{fig:pleiades_cmd_label}
\end{figure}

Arguably one of the simplest things to do with a cluster is to perform a multi-band photometric survey and create a colour-magnitude diagram (CMD) by plotting the magnitude versus the difference in magnitude in two different bandpasses. Fig.~\ref{fig:pleiades_cmd_label} shows the $V, V$--$I_{\rm{c}}$ CMD of the Pleiades with data taken from \cite{Stauffer07} in which several regions have been labelled as well as the global parameters one can estimate from stars in these regions. The Pleiades is arguably the best-studied (young) cluster with a rich sequence running from mid B-type stars down to the stellar/substellar boundary and so provides the ideal illustrative tool for the purposes of discussing the derivation of fundamental stellar parameters from clusters.

\subsection{Zero-age main-sequence}
\label{zams_distance}

At the age of the Pleiades the majority of stars have already settled onto the zero-age main-sequence (ZAMS) where they are steadily burning hydrogen in their cores. The ZAMS represents a temporally stable regime in which the stars essentially remain stationary with age in the CMD, and so whilst this region is not suitable for age estimates, it can however provide an estimate of the distance to the cluster; a prerequisite if one is interested in the bolometric luminosity of cluster members. So-called main-sequence fitting distances can be determined by comparing the sequence in question with either an empirical main-sequence relation based on photometry of nearby stars with known distances or via the use of theoretical stellar evolutionary models. There are of course issues related to both methods. For instance, empirical relations are typically mean relations computed from stars with a range of ages (of up to several Gyr) and so may not accurately reflect the positions of ZAMS stars in the CMD. The empirical sequence can be significantly brighter due to evolutionary effects and so will affect the apparent distance to the cluster (see e.g. \citealp{Littlefair10}). Similarly, stellar evolutionary models are not only susceptible to uncertainties in the underlying physics included, but must also be transformed from the theoretical Hertzsprung-Russell (H-R) plane to the observational CMD plane which requires the use of both a colour-$T_{\rm{eff}}$ relation to set the abscissa and bolometric correction (BC)-$T_{\rm{eff}}$ relations to set the ordinate. Regardless which method one adopts to infer the distance to a given cluster, the potential effects of both interstellar extinction and compositional differences must first be accounted for as both will act to modify the positions of stars in the CMD relative to unreddened and/or lower/higher metallicity stars.

Although it may seem that main-sequence fitting distances will be consigned to history in the age of \emph{Gaia}, it is worth remembering that there remains an outstanding discrepancy with regards to the \emph{Hipparcos} distance to the Pleiades and those estimated via main-sequence fitting (e.g. \citealp{An07}; cf. 120\,pc with 135\,pc) as well as other complementary methods (see e.g. \citealp{Soderblom05,Melis14}). Hence, although \emph{Gaia} will likely resolve this discrepancy, main-sequence fitting distances will continue to provide a simple, yet effective sanity check with regards to estimating cluster distances.

\subsection{Upper main-sequence}
\label{upper_ms_age_extinction}

The upper main-sequence (upper-MS) represents a region of the CMD which is not temporally stable i.e. as a function of age the evolution of stars becomes noticeable and can therefore be used to provide an age determination. When a star reaches the ZAMS it begins fusing hydrogen into helium in its core. Over time, the helium content of the core will increase and this leads to the star moving both redward and brighter in the CMD. Although subtle, this progression of the sequence between the ZAMS and terminal-age main-sequence (TAMS; where core hydrogen fusion ends) is a robust age indicator and the combination of stellar evolutionary models and sophisticated fitting techniques can provide reasonably well-constrained ages with statistically meaningful uncertainties (see e.g. \citealp{Naylor09, Bell13}).

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{Figs/upper.pdf}
 \caption{The effects of including rotation in stellar evolutionary models in the $M_{V}, B$--$V$ CMD. The isochrones shown are from the recent Geneva models by \cite{Ekstrom12}.}
   \label{fig:upper}
\end{figure}

The use of stellar evolutionary models means that any age determination will naturally be model-dependent and as such is prone to uncertainties inherent in the models. In the high-mass regime the main sources of uncertainty are the effects of stellar rotation, the degree of convective core overshoot and, to a lesser extent, the rate of mass-loss in early-type stars. The inclusion of rotation and/or convective core overshooting will act to increase the main-sequence lifetime of a star as a result of increased levels of hydrogen supplied to the core (see e.g. \citealp{Meynet00}). Fig.~\ref{fig:upper} illustrates the effects of including stellar rotation in evolutionary models. The models used in this demonstration are the recent Geneva models by \cite{Ekstrom12} for which the authors assume a fixed rotation rate of 40\% of break-up. It is evident from the $M_{V}, B$--$V$ CMD that the inclusion of rotation (the only difference between the two sets of isochrones shown in Fig.~\ref{fig:upper}) makes a significant difference to the position of the isochrones in the CMD, and ultimately the derived age (see e.g. \citealp{Brandt15b}). For example, ages determined from the main-sequence turn-off would be affected at the $\sim$30\% level. On the other hand, the transition from the ZAMS to TAMS is much less affected -- and of course there are a greater number of stars in this region compared to the post-main-sequence due to the shape of the mass function -- and so age determinations from this region are affected at only the $\sim$10\% level. The main issue with main-sequence ages of course is the relative paucity of stars (especially in younger clusters) and so whilst the resultant age may have a moderately small systematic uncertainty, the statistical uncertainty can be significant.

In addition to age estimates, the higher mass stars (late B-type and earlier) can also be used to calculate the extinction towards a given cluster. This is typically performed using multi-band photometry blueward of (and including) the $V$-band in either a colour-colour diagram or through the use of various reddening vectors to de-redden individual sources onto a given sequence (e.g. the Q-method; \citealp{Johnson53}). This can trivially and quickly allow one to ascertain whether the extinction is uniform across the cluster or whether it is spatially variable. Note that not all young clusters/star-forming regions contain high-mass stars (see e.g. Taurus), and so in such cases both photometry and spectroscopy of individual members is necessary i.e. the combination of an observed colour and a spectral type will permit one to determine the extinction of a given star.

\subsection{Pre-main-sequence}
\label{ms_turn-on_age}

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{Figs/pleiades_cmd_iso.pdf}
 \caption{The $V, V$--$I_{\rm{c}}$ CMD of the Pleiades with several sets of commonly used pre-MS isochrones overlaid. The upper reference refers to the stellar interior models, whereas the lower reference corresponds to the particular colour-$T_{\rm{eff}}$ and BC-$T_{\rm{eff}}$ relations used to transform the model into the observational plane.}
   \label{fig:pre-ms}
\end{figure}

Stars in the pre-main-sequence (pre-MS) phase are still contracting under the influence of gravity, and as such become noticeably fainter with age and hence provide an additional age diagnostic for a cluster. As with ages determined from the upper-MS, pre-MS ages are typically inferred via the use of stellar evolutionary models and are thus model-dependent. Unlike the upper-MS regime, however, there is a high degree of model dependency with pre-MS ages. Fig.~\ref{fig:pre-ms} illustrates this by showing several sets of publicly available pre-MS model isochrones which are commonly adopted in the literature overlaid on the $V, V$--$I_{\rm{c}}$ CMD of the Pleiades. The cluster parameters listed in the upper right of Fig.~\ref{fig:pre-ms} are all well-constrained and have been determined independently of fitting model isochrones in the CMD; namely the lithium depletion boundary (LDB) age from \cite{Barrado04b}, the VLBI distance of \cite{Melis14} and a reddening based on the mean extinction of \cite{Stauffer98a}. Fig.~\ref{fig:pre-ms} clearly shows a discrepancy between the observed Pleiades sequence in the CMD and the predicted colours/magnitudes of the models (having been transformed into the observational plane). Not only do none of the models match the sequence at cool temperatures (see e.g. \citealp{Bell12}), but furthermore they do so in a non-systematic way i.e. the derived age depends on which part of the sequence is fitted. Given the importance we place on the use of pre-MS model isochrones for estimating stellar ages, and by extension, timescales for important astrophysical phenomena in young (typically) less well-characterised clusters/associations, this is particularly perturbing.

If we are to use such models to infer ages from young star clusters it is apparent that some form of empirical correction to the BC-$T_{\rm{eff}}$ relation is required so as to fit an observed sequence in the CMD at a given age. Furthermore, if one wishes to use such ``semi-empirical'' models at significantly younger ages then any correction should include some form of surface gravity (log$\,g$) dependence (see e.g. \citealp{Bell14}). Note that even if such a correction is applied, the resultant ages from pre-MS models will still be heavily model-dependent and can differ by factors of 2--3 due to differences in the underlying input physics and parameters adopted in the models.

\subsection{Lithium depletion boundary}
\label{lithium_depletion_boundary}

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{Figs/ldb.pdf}
 \caption{The LDB in NGC~2547 defined by the clear discontinuity in the sequence between Li-poor and Li-rich members. The horizontal dashed lines denote constant luminosity loci with ages of 30, 40 and 50\,Myr. Figure from \protect\cite{Soderblom14}.}
   \label{fig:ldb}
\end{figure}

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{Figs/jackson_binks_ldb.pdf}
 \caption{The effects of including magnetic field-related phenomena on LDB ages. \emph{Left:} The effects of different star spot filling factors. Figure from \protect\cite{Jackson14b}. \emph{Right:} The effects of inhibiting convective flows. Figure from \protect\cite{Binks16}.}
   \label{fig:jackson_binks_ldb}
\end{figure*}

The basic premise underlying the LDB technique is that low-mass pre-MS stars (\ltilde\,0.4\,\msun) remain fully convective until they reach the ZAMS. As they contract, the core temperature increases until it is capable of burning Li at temperature of $\sim$\,3$\times$10$^{6}$\,K. Given the fully convective nature of the stars, Li-depleted material at the centre is brought up to the surface and Li-rich material from the surface brought to the centre where it becomes depleted, and complete Li depletion occurs in a fraction of a Kelvin-Helmholtz timescale. Fig.~\ref{fig:ldb} illustrates this rapid phase of Li destruction and the resultant sharp discontinuity in the $I, R$--$I$ CMD of NGC~2547 between stars which have contracted sufficiently to burn Li and those at slightly lower masses which have not. The LDB technique is evidently more telescope time intensive in the sense that one needs to first isolate potential members based on their positions in the CMD and then follow-up spectroscopy is required to identify the resonant Li feature at 6708\,\AA. The overriding benefit of this technique, however, is the high-degree of model-insensitivity. Unlike the ages derived from model isochrones in the CMD in which the various models simply do not agree, the different models agree remarkably well on the luminosity at which Li depletion occurs with systematic uncertainties at the $<$\,10\% level (see e.g. \citealp{Tognelli15}). Unfortunately, the LDB technique is practically only applicable in clusters with ages between 20 and 200\,Myr. At younger ages model dependency becomes increasingly problematic, reaching levels of $\sim$\,30\%, whereas at older ages the combined age and distance of potential target clusters means that the boundary itself simply becomes too faint. To date there are 10 clusters with LDB ages in the literature and these represent the necessary fiducial points against which other age-dating techniques can be validated in the same cluster (see e.g. \citealp{Mamajek14}).

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{Figs/hillenbrand_stassun_masses.pdf}
 \caption{\emph{Left:} Percentage difference between dynamically-determined stellar masses and those predicted by stellar evolutionary models for pre-MS stars in the H-R diagram. Figure from \protect\cite{Hillenbrand04}. \emph{Right:} Same as left-hand panel but showing only young eclipsing binaries. Figure from \protect\cite{Stassun14}.}
   \label{fig:hillenbrand_stassun_masses}
\end{figure*}

\subsection{The effects of magnetic fields}
\label{effects_magnetic_fields}

Although there is a very strong consensus amongst the various stellar evolutionary models as regards the luminosity of the onset of Li depletion, such models are well-known to be over simplistic representations of actual stars and do not include certain phenomena which we know to exist/occur in low-mass stars. Over the past couple of years several investigations have looked into the effects of including magnetic field-related phenomena in evolutionary models; namely the presence of star spots in the stellar photosphere and the inhibition of convective flows due to strong large-scale magnetic fields (see e.g. \citealp{Feiden14,Jackson14b,Somers15}). Both of these act to block outward flux from the star and hence slow the rate of contraction resulting in magnetic stars being cooler and larger at a given age than non-magnetic counterparts. As magnetic stars are cooler, it therefore takes longer to reach the necessary temperature in the core to burn Li and so the Li depletion timescale is extended. Fig.~\ref{fig:jackson_binks_ldb} illustrates the effects of including each of the aforementioned phenomena in stellar evolutionary models. Assuming a typical star spot coverage of $\sim$30\% (shown in the left-hand panel of Fig.~\ref{fig:jackson_binks_ldb}) or magnetic field strengths of a few kG (shown in the right-hand panel of Fig.~\ref{fig:jackson_binks_ldb}), the LDB age can increase by up to $\sim$20\% thereby suggesting that our current estimates represent lower limits to the cluster ages. Note the agreement between the two non-magnetic stellar evolutionary models in the right-hand panel of Fig.~\ref{fig:jackson_binks_ldb} which have been computed using different input physics and assumptions, that reiterates the statement concerning the high-level of model independence discussed in Section~\ref{lithium_depletion_boundary}.

\section{Young binary stars}
\label{binary_stars}

Binary stars represent a subset of stars which are sometimes referred to as ``benchmark'' or ``touchstone'' stars, in the sense that they permit tests of stellar evolutionary models at the most fundamental level by providing direct measurements of masses, radii, $T_{\rm{eff}}$ and bolometric luminosities. Of these properties, perhaps mass is of the greatest interest given that this is the primary input for stellar evolutionary models and essentially dictates the subsequent evolution of the star. Over the past decade several studies have collated the available dynamical information from young binary systems to test stellar evolutionary models, and in each case the same conclusion was reached; namely that modern evolutionary models are unable to accurately predict the properties of young, low-mass stars in binary systems (see e.g. \citealp{Hillenbrand04,Mathieu07,Stassun14}). The left-hand panel of Fig.~\ref{fig:hillenbrand_stassun_masses} shows that masses determined from positions in the H-R diagram are systematically underestimated for a given stellar $T_{\rm{eff}}$ and bolometric luminosity (i.e. radius). The study of \cite{Hillenbrand04}, however, was predominantly based on astrometric binaries with few eclipsing binaries in the sample. Over the past few years several such systems have been identified which led to a more recent study by \cite{Stassun14} who compared a newer generation of stellar evolutionary models to a selection of young eclipsing binaries only. Despite the inclusion of more up-to-date physics (e.g. opacities and molecular line lists), the models essentially fared no better than those investigated by \cite{Hillenbrand04}. For masses above 1\,\msun\ all models predict masses to within 10\% of the dynamical measurements, however below 1\,\msun\ the situation is much worse with mass errors of 50--100\%. Note a qualitative difference between the two panels in Fig.~\ref{fig:hillenbrand_stassun_masses}, namely that the \cite{Hillenbrand04} study found that the models tend to underestimate the mass, whereas the \cite{Stassun14} investigation found that the models overestimate the masses (with a potential dichotomy between the predicted primary and secondary masses). There are only three eclipsing binaries and one set of evolutionary models in common between the two studies for which very similar results were found.

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\textwidth]{Figs/kraus_binary.pdf}
 \caption{\emph{Left:} The primary and secondary components of UScoCTIO 5 in the H-R diagram. The isochrones and mass tracks are from the Dartmouth stellar evolutionary models of \protect\cite{Dotter08}. The blue lines denote the expected age of the system and the dynamically-determined mass of both components. \emph{Right:} Both components in the mass-radius diagram. The blue lines represent the expected age of the system and the measured luminosity of both components. Figure from \protect\cite{Kraus15}.}
   \label{fig:kraus_binary}
\end{figure*}

Recently, \cite{Kraus15} identified UScoCTIO 5, a known spectroscopic binary with spectral type M4.5, as an eclipsing binary with both eclipses apparent in K2 light curves. Fig.~\ref{fig:kraus_binary} shows the primary and secondary components of UScoCTIO 5 (which are almost of identical mass) in both the H-R and mass-radius diagrams. The blue lines in both panels represent the expected age, dynamical mass and measured luminosity of both components. It is clear from Fig.~\ref{fig:kraus_binary} that despite the on-going uncertainty with regards to the age one should adopt for particular stars in Upper Scorpius (see \citealp{Pecaut16}) the stellar evolutionary models (in 4 out of 5 cases) underestimate the masses. The notable exception to this are the recent PARSEC models of \cite{Chen14} which include an empirical correction to the outer boundary condition based on observations of low-mass stars in significantly older clusters and which appears to be an over-correction in the case of UScoCTIO 5.

There is of course a fundamental issue with using the H-R diagram to test stellar evolutionary models, which simply arises from the difficulty of measuring both the $T_{\rm{eff}}$ and bolometric luminosity. The former requires the use of model atmospheres which may still suffer from uncertainties regarding the underlying physics (especially in terms of missing sources of opacity) and the latter necessitates accurate and precise distance estimates (although upcoming \emph{Gaia} data releases will certainly help in this area). Thus, an even more fundamental test of stellar evolutionary models is to use the directly measured masses and radii and perform the comparison in the mass-radius diagram. The right-hand panel of Fig.~\ref{fig:kraus_binary} demonstrates that the models predict radii which are too small for the known luminosity and mass of the stars, which subsequently implies that the model $T_{\rm{eff}}$ are too high. This $T_{\rm{eff}}$ offset could be due to a number of potential issues including underlying issues with the prescription of convection in low-mass stars, missing opacities as well as a miscalibration of the spectral type-$T_{\rm{eff}}$ scale.

\begin{figure*}[t]
\centering
\includegraphics[width=0.7\textwidth]{Figs/binary.pdf}
 \caption{The mass-radius diagram of young eclipsing binaries in which both the primary and secondary components have masses of \ltilde\,1.4\,\msun\ with the recent stellar evolutionary models of \cite{Baraffe15} overlaid.}
   \label{fig:binary}
\end{figure*}

Fig.~\ref{fig:binary} shows the mass-radius diagram for eclipsing binaries with ages of the Pleiades and younger in which both the primary and secondary components have masses of \ltilde\,1.4\,\msun, and which are compared against the recent stellar evolutionary models of \cite{Baraffe15}. There are a couple of points worth noting. First, although the number of systems has increased over the past few years, there is still a paucity of young eclipsing binaries with which to rigorously test evolutionary models. It is likely that the continuing K2 mission will help identify more such systems in the coming years. For example, if one can identify a number of young eclipsing binaries spanning a range of masses at a given age, across a range of ages, this will then not only permit stringent tests of the models, but may also perhaps allow us to constrain the uncertain physics inherent to them. Second, there appears to be a mixture of binary systems which are coeval and systems which are not. Such instances have previously been discussed in the literature (see e.g. \citealp{Gennaro12}), however it is worth noting that even the comparison of data and models in the mass-radius diagram is not strictly a fair comparison. For example, one is comparing low-mass, likely magnetically active stars in a binary system (with an increased potential for strong dynamical interactions at the earliest phases) with standard single-star, non-magnetic models. This could further be complicated by possible episodic accretion histories which may have dramatic effects on the internal structure of the stars (see e.g. \citealp{Audard14}).

\section{Spectroscopic surveys}
\label{spectroscopic_surveys}

\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{Figs/jackson_radii.pdf}
 \caption{\emph{Top:} Stellar radius versus bolometric luminosity for stars in three young clusters; NGC~2516 (140\,Myr), NGC~2547 (35\,Myr) and NGC~2264 (5\,Myr). The black squares denote the averaged radii in specific luminosity bins. The following stellar evolutionary models are shown for comparison: BHAC15 (black solid; \citealp{Baraffe15}), Dartmouth (blue dotted; \citealp{Dotter08}), Dartmouth modified to include magnetic fields (blue dashed; \citealp{Feiden15}), YREC (red dot-dashed; \citealp{Demarque08}), and YREC modified with an effective star spot coverage of 30\% (red solid; \citealp{Somers15}). The green shaded region represents the expected luminosities over which stars develop radiative cores according to the BHAC15 models; stars to the left have radiative cores, whereas those to the right are fully convective. \emph{Bottom:} The percentage difference between the measured radii and those predicted by the BHAC15 models. The blue dashed and red solid lines represent the percentage difference between the non-magnetic and modified Dartmouth and YREC evolutionary models respectively. Figure from \protect\cite{Jackson16}.}
   \label{fig:jackson_radii}
\end{figure*}

In the last few years the \emph{Gaia}-ESO and APOGEE/IN-SYNC surveys have collected large-scale spectroscopic samples of stars in young clusters. The primary benefit of such surveys is the homogeneous datasets which have been produced for the community which permit the investigation of common features and peculiarities of different clusters in a self-consistent and standardised way. This represents a marked shift from bringing together the hodgepodge of smaller surveys by different groups looking at different subsets of stars in a cluster and which have not only been collected using a variety of instruments, but have also been reduced and analysed in a heterogeneous way.

\begin{figure*}[hbtp!]
\centering
\includegraphics[width=\textwidth]{Figs/cottaar_ic348.pdf}
 \caption{Distributions of three (almost) independent stellar radius estimates for members of IC~348. The blue points and histograms show the observed distributions. The black distributions represent the expected distributions accounting for observational effects which could introduce an observed spread in stellar radii including binarity, measurement uncertainties and projection effects. The coloured distributions are the same as the black distributions but also include a Gaussian intrinsic stellar radius spread with a width of 15\% (orange), 25\% (red) and 35\% (magenta) around the median stellar radius at the corresponding $T_{\rm{eff}}$. Figure from \protect\cite{Cottaar14}.}
   \label{fig:cottaar_ic348}
\end{figure*}

\subsection{Inflated radii and radii spreads}
\label{inflated_radii_radii_spreads}

Observations have long suggested that the radii of short-period eclipsing binaries are inflated with respect to those predicted by standard stellar evolutionary models (see e.g. \citealp{Kraus11}). Combining rotational periods with projected equatorial velocities from the \emph{Gaia}-ESO survey, \cite{Jackson16} estimated the average radii of cluster members as a function of luminosity and age to investigate whether similarly inflated radii are observed in young, low-mass, rapidly-rotating stars. Fig.~\ref{fig:jackson_radii} shows the average radii of cluster members as a function of luminosity for three young clusters; NGC~2516 (140\,Myr), NGC~2547 (35\,Myr) and NGC~2264 (5\,Myr), which are compared against several sets of stellar evolutionary models (including modified versions of both the Dartmouth and YREC models which incorporate a prescription for magnetic fields and a star spot coverage of 30\% respectively, see \citealp{Feiden15,Somers15}). The lower panels of Fig.~\ref{fig:jackson_radii} display the percentage difference between the observed and predicted radii from the BHAC15 models \citep{Baraffe15}, whereas the blue dashed and red solid lines represent the percentage difference between the non-magnetic and modified Dartmouth and YREC evolutionary models respectively. In short, the radii of young cluster members are larger than those predicted by standard stellar evolutionary models and is more pronounced in fully convective low-mass stars in which it can reach levels of $\sim$30\%. With the modified evolutionary models currently available, it appears as though strong surface magnetic fields (exceeding 2.5\,kG), star spots covering $\sim$30\% of the photosphere or, more likely, a combination of both is necessary to explain the observed data.

The \emph{Gaia}-ESO and IN-SYNC spectroscopic surveys are providing high-precision fundamental stellar parameters (e.g. $\sigma_{\rm{RV}}$\,\ltilde\,0.2\,km\,s$^{-1}$, $\sigma_{\mathrm{log}\,g}$\,\less\,0.1\,dex, $\sigma_{T_{\rm{eff}}}$\,\less\,100\,K) and such precision is allowing us to investigate the constituent members of young star clusters in unprecedented detail. The first cluster observed within the framework of the APOGEE/IN-SYNC ancillary programme was IC~348 (see \citealp{Cottaar14}). In brief, the IN-SYNC programme has acquired thousands of high-resolution $H$-band spectra for thousands of pre-MS stars in young clusters. From these spectra they extract various parameters including $T_{\rm{eff}}$, surface gravity (log\,$g$) and rotational velocities ($v$ sin\,$i$) via comparison with stellar atmospheric models. Fig.~\ref{fig:cottaar_ic348} shows the observed distributions for three (almost) independent stellar radius estimates for members of IC~348; namely the projected stellar radius, surface gravity and extinction-corrected $J$-band magnitude. From Fig.~\ref{fig:cottaar_ic348} it is clear that compared to the expected distributions, which account for observational effects that could introduce an observed radius spread including binarity, measurement uncertainties and projection effects, the observed distributions are significantly wider. \cite{Cottaar14} find a best-fit to the observed distributions by combining this expected distribution with a Gaussian intrinsic stellar radius spread of width 25\% around the median stellar radius at the corresponding $T_{\rm{eff}}$. Furthermore, \citeauthor{Cottaar14} argue that all three radius diagnostics are correlated with the brighter stars tending to have significantly lower surface gravities and larger projected stellar radii, and that uncertainties in the derived parameters are not responsible for the observed spread.

The simplest explanation for an intrinsic spread in stellar radii in IC~348 is that this is a result of an intrinsic age spread within the cluster i.e. the younger stars have not contracted as much as the older stars in the same cluster and hence have correspondingly lower surface gravities and larger projected stellar radii. The exact age spread, however, is dependent upon the mean age of the IC~348 which lies between 3 and 6\,Myr. The main uncertainty in the age of the cluster is the uncertainty on its distance (ranging from 220 to 350\,pc; see e.g. \citealp{Herbst08}) and so hopefully this is another instance in which \emph{Gaia} will provide relief in the coming years. Assuming the older age of 6\,Myr, this corresponds to an upper limit on the age spread of 8\,Myr. Note, however, that \cite{Cottaar14} also find that the more rapid rotators have larger stellar radii and thus it is likely that the spread in stellar radii is not solely due to an intrinsic age spread, but may also include contributions from different accretion histories and/or different levels of magnetic activity amongst cluster members. The ability to disentangle these potential effects from one another to infer any genuine intrinsic age spread is currently beyond our means, however if intrinsic age spreads of a few Myr are common amongst the youngest clusters/star-forming regions (see also \citealp{DaRio16}) then this could have serious implications, especially regarding their use as fiducial age ``points'' in understanding, for example, circumstellar disc dissipation timescales, and hence planet formation timescales, as by definition the cluster does not have a single age, but a range.

\subsection{Kinematic substructure}
\label{kinematic_substructure}

The low-mass stellar population associated with the Wolf-Rayet binary $\gamma^{2}$ Velorum was the first young ``cluster'' target observed as part of the \emph{Gaia}-ESO survey (see \citealp{Jeffries14}). Fig.~\ref{fig:jeffries_gamma_vel} shows the binned radial velocity histogram for $\gamma^{2}$ Velorum members. It is apparent that the best-fit model consisting of a single Gaussian component, which includes a fraction of unresolved binaries according to the results of the recent \citealp{Raghavan10} survey, represents a poor fit to the data. By contrast, the best-fit model comprised of two Gaussian components, each with an unresolved binary population, is a much better fit. The small uncertainties on the measured radial velocities are such that \citeauthor{Jeffries14} were able to unambiguously identify two distinct kinematic populations (hereafter referred to as Populations A and B, each of which contains approximately equal numbers of stars); one with a very narrow intrinsic width and the other with a much broader width (cf. $\sigma_{\rm{RV, A}}$\,\equal\,\,0.34\,\plusminus\,0.16 and $\sigma_{\rm{RV, B}}$\,\equal\,1.60\,\plusminus\,0.37\,km\,s$^{-1}$) and which is offset from the first by 2\,km\,s$^{-1}$. \citeauthor{Jeffries14} also identified several other differences between the two Populations, notably: i) based on levels of Li depletion, Population~A is about 1--2\,Myr older than Population~B and ii) Population~A appears to be in viral equilibrium and likely represents the bound remnant of an initially larger cluster which formed in the denser region of the Vela~OB2 association, whereas Population~B comprises a dispersed population of unbound stars which probably formed in the less dense regions of the association. Interestingly, the low-mass stars of Populations~A and B appear to be several Myr older than $\gamma^{2}$ Velorum, thus suggesting a scenario in which it formed after the bulk of the low-mass population and which possibly resulted in the termination of star formation in the region via the expulsion of gas and dust.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{Figs/jeffries_gamma_vel_orig.pdf}
 \caption{Binned radial velocity histogram for members of the $\gamma^{2}$ Velorum cluster. \emph{Top:} Best-fit for a model consisting of a single Gaussian component including unresolved binaries. \emph{Bottom:} Same as the upper panel but instead the model consists of two Gaussian components each with an unresolved binary population. Figure from \protect\cite{Jeffries14}.}
   \label{fig:jeffries_gamma_vel}
\end{figure}

\section{Summary of conclusions}
\label{concluding_thoughts}

Below I briefly reiterate the main conclusions from this contribution.

\begin{enumerate}
\item The colour-magnitude diagram of a given cluster can provide several global parameters shared by constituent members (including age, distance and the presence/uniformity of interstellar extinction), although one should be aware of the underlying uncertainties as regards the use of stellar evolutionary models and carefully assess the pedigree of adopted empirical relations.
\item Model-dependent estimates of low-mass pre-MS stellar parameters are unreliable. Specifically, below 1\,\msun\ the models tend to underestimate the mass of a given star based on its position in the H-R diagram as well as predict radii which are too small based on its position in the mass-radius diagram.
\item There is tentative evidence that the introduction of magnetic field-related phenomena (such as star spots and/or the inhibition of convective flows) may help to resolve the discrepancy between dynamically-determined parameters and those predicted by stellar evolutionary models.
\item Recent spectroscopic survey results have demonstrated that clusters are more complex entities than previously thought (e.g. age spreads and kinematic substructure) and the continuing \emph{Gaia}-ESO and APOGEE/IN-SYNC surveys will only highlight further examples of this and hence continue to shape our understanding of the formation and early evolutionary stages of young clusters.
\end{enumerate}

\section*{Acknowledgments}
{I would like to thank the organisers of the splinter session for inviting me to present this review as well as the SOC/LOC who together provided an excellent scientific programme and an exceptionally ran conference.}

\bibliographystyle{mn3e3}


\section{Introduction}
\label{sec:introduction}

Two major directions are being explored to enhance radio frequency (RF) spectrum access, namely, reusing currently under-utilized frequency bands below roughly 6~GHz through dynamic spectrum access (DSA) and exploiting spectrum opportunities in higher frequencies, e.g., millimeter-wave (mmWave) bands \cite{mingmingcaiJNLAsilomar2015, el2014spatially}. The total amount of spectrum that can be shared below 6~GHz is limited to hundreds of MHz, whereas mmWave bands can provide spectrum opportunities on the order of tens of GHz.  However, the signals at mmWave frequencies experience significantly higher path loss than signals below 6~GHz \cite{rappaporttheodores2002}.

Beamforming with a phased array having a large number of antennas is a promising approach to compensate for higher attenuation at mmWave frequencies while supporting mobile applications \cite{huang2015joint,hur2013millimeter, heath2015overview}.
First, the short wavelength of mmWave allows for integration of a large number of antennas into a phased array of small size suitable for a mobile wireless device.
Second, switched beamforming applies different sets of parameters to the phased array, allowing the beam direction to be electronically and quickly controlled, which enables steering for different orientations and tracking for mobility.  In this paper, we consider switched analog beamforming with one RF chain using a phased array with phase shifters.  For hybrid architectures, phased arrays can also be implemented by either phase shifters or switches \cite{mendez2015hybrid}.

Traditionally in a communication system with analog beamforming, the sets of phase shifter values in a phased array are designed at a specific frequency, usually the carrier frequency, but applied to all frequencies within the transmission bandwidth due to practical hardware constraint. Phase shifters are relatively good approximations to the ideal time shifters for narrowband transmission; however, this approximation breaks down for wideband transmission if the angle of arrival (AoA) or angle of departure (AoD) is away from broadside, because the required phase shifts are frequency-dependent.
The net result is that beams for frequencies other than the carrier ``squint'' as a function of frequency \cite{mailloux2005phased} in a wide signal bandwidth. This phenomenon is called beam squint in a wideband system \cite{SeyedGarakkoui2011BeamSquinting}. As we will see in the context of orthogonal division multiplexing (OFDM) modulation with multiple subcarriers, beam squint translates into array gains for a given AoA or AoD that vary with subcarrier index.  To eliminate or reduce beam squint with hardware approaches, true time delay and phase improvement schemes have been proposed  \cite{longbrake2012true, liu2013minimize}; however, these approaches are unappealing for mobile wireless communication due to combinations of high implementation cost, significant insertion loss, large size, or excessive power consumption.

Covering a certain range of AoA or AoD in switched beamforming requires a ``codebook'' consisting of multiple beams, with each beam determined by a set of beamforming phases \cite{song2015codebook, cui2014evolution}. Usually, a minimum array gain is required for the codebook to cover the entire desired angle range \cite{hehierarchical}, for example, 3~dB below the array's maximum gain. A natural extension for a wideband system is to design the beamforming codebook to guarantee the minimum gain for all frequencies, even in the presence of beam squint.  To the best of our knowledge, only \cite{liu2013minimize} proposes a phase improvement scheme to reduce the effect of beam squint in beamforming codebook, but with high implementation cost of additional phase shifters and bandpass filters in mmWave frequency.
By contrast, we propose developing denser codebooks to compensate for beam squint at the system level.

The main contributions of this paper are twofold.  First, we model the beam squint for a uniform linear array (ULA) in a wireless communication system, and show that the usable beamwidth for each beam effectively decreases because of the beam squint. Therefore, the number of beams in the codebook must increase to counter this effect.
Second, we develop a codebook design algorithm to enforce minimum array gain for all subcarriers in the wideband system. The bandwidth of the system is limited for a fixed number of antennas in the phased array; conversely, the number of antennas in the phased array is limited for a fixed bandwidth. The algorithm and conclusions also apply to frequencies other than mmWave and acoustic communication.

\section{System Model}
\begin{figure}
\begin{center}
\includegraphics[width=85 mm]{figs/array-structure}
\end{center}
\caption{Structure of a uniform linear array (ULA) with analog beamforming using phase shifters.  The distance between the adjacent antennas is $d$, and $\theta$ denotes either the angle-of-arrival (AoA) for reception or the angle-of-departure (AoD) for transmission.}
\label{fig:array-structure}
\end{figure}

We consider a ULA as shown in Fig.~\ref{fig:array-structure}. It could be the phased array of either a transmitter or receiver. In the array, there are $N$ identical and isotropic antennas, each connected to a phase shifter. The distances between two adjacent antenna elements are the same and are all denoted by $d$. For simplicity, we also assume the phases of the phase shifters are continuous without quantization. The AoD or AoA $\theta$ is the angle of the signal relative to the broadside of the array, increasing counterclockwise.
From \cite{orfanidis2002electromagnetic}, the ULA response vector is
\begin{align}
{\mathbf{a}}\left( \theta \right) =
\left[ 1,{e^{j2\pi {\lambda ^{ - 1}}d\sin \theta}}, \ldots ,{e^{j2\pi {\lambda ^{ - 1}}nd\sin\theta }}, \ldots , \right. \nonumber \\
 \left. {e^{j2\pi {\lambda ^{ - 1}}\left( {N - 1} \right)d\sin \theta}} \right]^T, \ \theta \in {\left[ -\frac{\pi}{2}, \frac{\pi}{2} \right]},
\label{eq:receive-array-response-linear-array}
\end{align}
where $\lambda$ is the signal wavelength, and superscript $\left( \cdot \right) ^T$ indicates transpose of a vector. The effect of the phase shifters can be modeled by the beamforming vector
\begin{align}
{\mathbf{w}} = \left[ {e^{j\beta _1},\dots,e^{j\beta_n},\dots,e^{j\beta_N}} \right]^T,
\label{eq:beam-forming-vector}
\end{align}
where $\beta_n$ is the phase of $n$th phase shifter connected to the $n$th antenna element in the array.

From \cite{balanis2016antenna}, the array factors or array gains for both the transmitter and receiver have the same expression
\begin{align}
g \left({\mathbf{w}},\theta \right)=\frac{1}{\sqrt {N}} {\mathbf{w}}^H  {\mathbf{a}}\left( \theta \right) = \frac{1}{\sqrt {N}} \sum\limits_{n = 1}^{N} {{e^{j\left[ {2\pi {\lambda ^{ - 1}}\left( {n - 1} \right)d\sin \theta - {\beta _n}} \right]}}}.
\label{eq:array-gain}
\end{align}
where superscript $\left( \cdot \right) ^H$ denotes Hermitian transpose.

\section{Modeling of Beam Squint in a Phased Array}

The analysis of the effect of beam squint is the same for a transmit or a receive array, so without loss of generality, we focus on the case of a receive array. Define the \emph{beam focus angle} $\theta_0\left(\mathbf{w}\right)$ to be the AoA for which the beam has the highest gain for a given beamforming vector $\mathbf{w}$. The beamforming vector $\mathbf{w}$ is traditionally designed for the carrier frequency. From (\ref{eq:array-gain}), to achieve the highest array gain for a beam focus angle $\theta_0$ among all possible $\mathbf{w}$, the phase shifters should be
\begin{align}
{\beta _n} \left( \theta_0 \right)=  2\pi {\lambda_c^{ - 1}}\left( {n - 1} \right)d\sin \theta_0, \ n=1,2,\dots,N,
\label{eq:phases-required-for-fine-beams}
\end{align}
where $\lambda_c$ is the wavelength of the carrier frequency. We call the beam with phase shifts described in (\ref{eq:phases-required-for-fine-beams}) a \emph{fine beam}, and we focus on the analysis of fine beams throughout the paper.

Based on (\ref{eq:array-gain}), the array gain for a signal with frequency $f$ and AoA $\theta$ using beam focus angle $\theta_0$ is
\begin{align}
g \left({\theta_0,  \theta, f}\right) = \frac{1}{{\sqrt {{N}} }}   \sum\limits_{n = 1}^N {{e^{j\left[ {2\pi {f c^{-1}}\left( {n - 1} \right)d\sin {\theta} - {\beta _n}\left(\theta_0 \right)} \right]}}},
\label{eq:ms-phase-combiner-implementation-freq}
\end{align}
where $c$ is the speed of light.

\begin{figure}
\begin{center}
\includegraphics[width=88 mm]{figs/Beam_plot_several_frequencies}
\end{center}
\caption{Example of beam patterns varying for different frequencies in a wideband system with beam focus angle $\theta_0=\pi/6$, $N=16$ antennas, antenna spacing $d=\lambda_c/2$, carrier frequency$f_c=73$~GHz, minimum frequency $f_{min}=65.7$~GHz, and maximum frequency $f_{max}=80.2$~GHz. The bandwidth is chosen to be quite large for purposes of illustration.}
\label{fig:beam-plot-several-frequencies}
\end{figure}

For a fixed array and beamforming vector, an example of the beam squint patterns for different frequencies is shown in Fig.~\ref{fig:beam-plot-several-frequencies}. The beam focus angle for the carrier frequency is $\pi/6$. The array gain for an AoA $\theta$ are different for different frequencies; in particular, the array gains of $f_{max}$ and $f_{min}$ for the beam focus angle $\theta_0$ are both smaller than that of $f_{c}$.

For a wideband system, we express a subcarrier's absolute frequency $f$ relative to the carrier frequency $f_c$ as
\begin{align}
f=\xi f_c,
\label{eq:freq-coefficient}
\end{align}
where $\xi$ is the ratio of the subcarrier frequency to the carrier frequency. Suppose the baseband bandwidth of the signal is $B$. Then $f \in \left[ f_c-\frac{B}{2}, f_c+\frac{B}{2} \right]$. Define the fractional bandwidth
\begin{align}
b:={B}/{f_c},
\label{eq:fractional-bandwidth}
\end{align}
with $b>0$, so that $\xi \in \left[1-b/2, 1+b/2\right]$. As a specific example, $\xi$ varies from 0.983 to 1.017 in a system with 2.5~GHz bandwidth at 73~GHz carrier frequency.

The expression in (\ref{eq:ms-phase-combiner-implementation-freq}) can be transformed  into
\begin{align}
g \left({\theta_0,  \theta_c, \xi}\right) &= \frac{1}{{\sqrt {{N}} }}  \sum\limits_{n = 1}^{N} {{e^{j\left[ {2\pi {\xi f_c c^{-1}}\left( {n - 1} \right)d\sin {\theta_c}- {\beta _n}\left(\theta_0 \right)} \right]}}} & \nonumber\\
&=\frac{1}{{\sqrt {{N}} }}  \sum\limits_{n = 1}^{N} {{e^{j\left[ {2\pi { \lambda_c ^{-1}}\left( {n - 1} \right)d \left(\xi \sin {\theta_c} \right) - {\beta _n}\left(\theta_0 \right)} \right]}}},
\label{eq:ms-phase-combiner-implementation-ratio}
\end{align}
where $\theta_c$ is AoA for the carrier frequency.

Let ${\theta^{\prime}}$ be the equivalent AoA for the subcarrier with coefficient $\xi$. Comparing (\ref{eq:ms-phase-combiner-implementation-ratio}) to (\ref{eq:ms-phase-combiner-implementation-freq}), $\theta^{\prime}$ can be expressed as
\begin{align}
\theta^{\prime} = \arcsin \left( {\xi \sin \theta_c} \right).
\end{align}
The variation of beam patterns for different frequencies in the wideband system can therefore be interpreted as the variation of AoAs in the same angle response for different subcarriers.

Usually, the beamforming array is designed such that $d=\lambda_c/2$. From (\ref{eq:phases-required-for-fine-beams}) and (\ref{eq:ms-phase-combiner-implementation-ratio}),
\begin{align}
g & \left({\theta_0,  \theta_c, \xi}\right) = \frac{1}{{\sqrt {{N}} }}  \sum\limits_{n = 1}^{N} {{e^{j\left[ {\pi \left( {n - 1} \right) \left(\xi \sin {\theta_c}-\sin {\theta_0} \right) } \right]}}} \nonumber \\
&= \frac{{\sin \left( {\frac{{N\pi }}{2}\left( {\xi \sin \theta_c  - \sin \theta_0 } \right)} \right)}}{{\sqrt N \sin \left( {\frac{\pi }{2}\left( {\xi \sin \theta_c  - \sin \theta_0 } \right)} \right)}} e^{j \frac{\left(N-1\right) \pi} {2} \left( {\xi \sin \theta_c- \sin \theta_0} \right)}.
\label{eq:ms-phase-combiner-implementation-ratio-3}
\end{align}
Define
\begin{align}
g \left(x \right) = \frac{{\sin \left( {\frac{{N\pi x}}{2}} \right)}}{{\sqrt N \sin \left( {\frac{\pi x }{2}} \right)}} e^{j \frac{\left(N-1\right) \pi x} {2}}.
\label{eq:ms-array-gain-one-parameter}
\end{align}
The array gain for frequency with $\xi$ in (\ref{eq:ms-phase-combiner-implementation-ratio-3}) becomes $g \left(\xi \sin \theta_c  - \sin \theta_0 \right)$. Fig.~\ref{fig:Beam_shape_plot_fine_beam_psi} illustrates an example plot of $\left|g \left(x \right) \right|$. The main lobe is the lobe with the highest gain.

\begin{figure}

\begin{center}
\includegraphics[width=88 mm]{figs/Beam_shape_plot_fine_beam_psi}
\end{center}
\caption{Example of $\left|g \left(x \right) \right|$ from (\ref{eq:ms-array-gain-one-parameter}) for a fine beam with $N=16$ antennas and antenna spacing $d=\lambda_c/2$.}
\label{fig:Beam_shape_plot_fine_beam_psi}
\end{figure}

Let
\begin{align}
\psi=\sin \theta
\end{align}
to convert the angle $\theta$ space to the $\psi$ space. Correspondingly, $\psi_c=\sin \theta_c$, $\psi_0=\sin \theta_0$, and $\psi_0 \left(\mathbf{w}\right)=\sin{\theta_0 \left(\mathbf{w}\right)}$. The array gain of a fine beam in (\ref{eq:ms-phase-combiner-implementation-ratio-3}) then becomes
\begin{align}
g \left(\xi \psi_c - \psi_0\right) =\frac{{\sin \left( {\frac{{N\pi }}{2}\left( {\xi  \psi_c  -  \psi_0 } \right)} \right)}}{{\sqrt N \sin \left( {\frac{\pi }{2}\left( {\xi \psi_c  - \psi_0 } \right)} \right)}} e^{j \frac{\left(N-1\right) \pi} {2} \left( {\xi \psi_c-  \psi_0} \right)},
\label{eq:ms-phase-combiner-implementation-ratio-4}
\end{align}
where $\varphi_c, \varphi_0 \in \left[-1,1\right]$. The maximum array gain is
\begin{align}
g_m=\mathop {\max }\limits_{\psi_c \in \left[-1,1\right]}   \left| g \left({\psi_c-\psi_0}\right) \right|=\sqrt{N}.
\label{maximum-array-gain}
\end{align}

For simplicity, we study the effect of beam squint in $\psi$ space. From (\ref{eq:ms-phase-combiner-implementation-ratio-4}), the effect of beam squint increases if the AoA is away from the beam focus angle, and it tends to increase as the beam focus angle increases.

\section{Codebook Design}
One widely used codebook design criterion is to require the array gain to be larger than a threshold $g_t$ for all subcarriers \cite{hehierarchical}. The threshold may vary for different system designs. In communication systems, the beam focus angle ${\psi _0}\left(\mathbf{w}\right)$ is determined by the beamforming vector $\mathbf{w}$. Define the beam coverage of $\mathbf{w}$ as the set
\begin{align}
\mathcal{R}\left( \mathbf{w} \right):
=&\left\{ \psi_c \in \left[ -\psi_m, \psi_m \right]: \right. \nonumber \\
&\left.\ \min_{\xi \in \left[1-b/2, 1+b/2\right] } \left| {g \left( { \xi \psi_c  -  {\psi _0}\left( {{{\mathbf{w}}}} \right) } \right)}\right| \ge {g_t} \right\},
\label{beam-range-definition}
\end{align}
where $0<\psi_m \le 1$ is the maximum targeted covered AoA/AoD of the codebook. The beam coverage of some beams is not contiguous in $\psi$ if $\psi_0$ is around -1 or 1. We only consider the beam whose beam range is contiguous.

Define the effective beamwidth $\Delta \psi$ as the measure of $\mathcal{R}\left(\mathbf{w}\right)$, i.e.,
\begin{align}
\Delta \psi  := \int_{ - {\psi _m}}^{{\psi _m}} {{I_{{\mathcal R}\left( {\mathbf{w}} \right)}}} \left( \psi  \right) d\psi,
\label{beam-width-definition}
\end{align}
where
\begin{align}
{I_{{\mathcal R}\left( {\mathbf{w}} \right)}}\left( \psi  \right) = \left\{
{\begin{array}{*{20}{l}}
   {1,\ \text{if}\ \psi  \in {\mathcal R}\left( {\mathbf{w}} \right)}  \\
   {0,\ \text{otherwise} }. \\
\end{array}} \right.
\end{align}

A single beam can only cover a small range of angles in space. In communication systems with beamforming, multiple beams with different beam focus angles are utilized to cover the target angle range $\left[ -\psi_m, \psi_m \right], 0<\psi_m \le 1$, forming a codebook $\mathcal{C}$. Suppose there are $M$ beams in the codebook to cover $\left[ -\psi_m, \psi_m \right], 0<\psi_m \le 1$. The codebook $\mathcal{C}$ has size $M=\left| \mathcal{C} \right|$ and can be defined as
\begin{align}
{\mathcal C} := \left\{ {{{\bf{w}}_m},m = 1,2,\dots,M:\bigcup\limits_{m = 1}^M {{\cal R}\left( {{{\mathbf{w}}_m}} \right)}  = \;\left[ { - {\psi _m},{\psi _m}} \right]} \right\}.
\label{codebook-definition}
\end{align}

Denote $\mathbf{\Omega}$ as the set of all codebooks that meet the requirements of (\ref{codebook-definition}). There is an infinite number of codebooks in $\mathbf{\Omega}$. Among them, we focus on the ones that achieve the minimum codebook size. Minimum codebook size is defined as
\begin{align}
{\left| \mathcal{C} \right| }_{min} =  \mathop {\min }\limits_{\mathcal{C} \in \mathbf{\Omega}} {{\left| \mathcal{C} \right| }}.
\end{align}

The threshold $g_t$ can be different for different system requirements. For illustration, we choose $g_t$ such that
\begin{align}
g_t={g_m}/ \sqrt{2},
\label{threshold-definition}
\end{align}
where $g_m$ is defined in (\ref{maximum-array-gain}). The analysis of other thresholds $g_t$ is similar.

\subsection{Codebook Design without Considering Beam Squint}
\begin{algorithm}[t]
\caption{Codebook Design without Beam Squint}
\label{codebook-design-no-beam-squint}
\begin{algorithmic}[1]
\State Calculate ${\left| \mathcal{C} \right| }_{min}$ from (\ref{minimum-codebook-size-no-beam-squint})
\If {${\left| \mathcal{C} \right| }_{min}$ is odd}
    \State Let $\psi_0=0$, $k=1$
    \State Align the first beam so that $\psi_0=0$
\Else
    \State Let $\psi_0= -\Delta\psi_{3\text{dB}}/2$, $k=0$
\EndIf
\State $\psi_0=\psi_0+\Delta\psi_{3\text{dB}}$
\While {$k < {\left| \mathcal{C} \right| }_{min}$}
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $\psi_0$
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $-\psi_0$
    \State $\psi_0=\psi_0+\Delta\psi_{3\text{dB}}$
\EndWhile
\end{algorithmic}
\end{algorithm}

Define beam edges as the points with $\psi$ in the main lobe whose array gain is $g_t$. There are two beam edges, one on the left side and one on the right side, denoted as $\psi_l$ and $\psi_r$, respectively. The selected $g_t$ in (\ref{threshold-definition}) is 3~dB below the maximum array gain $g_m$. The corresponding beamwidth is called the 3~dB beamwidth, denoted as $\Delta\psi_{3\text{dB}}$ in $\psi$ space. For consecutive beam coverage, the beamwidth is the distance between the two beam edges, i.e.,
\begin{align}
\Delta\psi_{3\text{dB}}=\psi_r-\psi_l.
\end{align}
According to \cite{orfanidis2002electromagnetic},
\begin{align}
\Delta\psi_{3\text{dB}}=\frac{1.772}{N},
\label{beam-width-N}
\end{align}
which is a constant independent of $\psi_0$. It is worth stressing that the beamwidth in $\theta$ space is not fixed, and depends on $\theta_0$.

To design a codebook with the minimum codebook size, we simply tile the beam coverages over $\left[ -\psi_m, \psi_m\right]$ without overlap.  If the effect of beam squint is not considered, all beams have the same beamwidth $\Delta\psi_{3\text{dB}}$ as shown in (\ref{beam-width-N}), and the minimum codebook size
\begin{align}
{\left| \mathcal{C} \right| }_{min}= \ceil[\bigg]{\frac{2\psi_m}{\Delta\psi_{3\text{dB}}}},
\label{minimum-codebook-size-no-beam-squint}
\end{align}
where $\ceil[\big]{\cdot}$ denotes the well-known ceiling function.

Even for the same ${\left| \mathcal{C} \right| }_{min}$, there could be multiple codebook designs. Among them, we are more interested in the ones that are symmetric relative to broadside. A codebook design method for minimum codebook size without considering the beam squint is shown in Algorithm~\ref{codebook-design-no-beam-squint}.

\subsection{Beam Pattern with Beam Squint}
As analyzed in the previous section, the array gains of some subcarriers are smaller than that of the carrier frequency when $\psi_c \ne 0$. To guarantee a minimum array gain $g_t$ for all the subcarriers, the beam range $\mathcal{R}\left( \mathbf{w} \right)$ is reduced.

Since $\mathbf{w}$ determines the beam focus angle $\psi_0$, $\mathcal{R}\left( \mathbf{w} \right)$ and $\mathcal{R}\left( \psi_0 \right)$ can be used interchangeably. Define
\begin{align}
&{\psi _{c,l}\left(\psi_0\right)}=\min \mathcal{R}\left( \psi_0 \right), \\
&{\psi _{c,r}\left(\psi_0\right)}=\max \mathcal{R}\left( \psi_0 \right).
\end{align}
First, we study the case for $\psi_c>0$, which gives us
\begin{align}
&{\psi _{c,l}\left(\psi_0\right)}\left( {1 - b/2} \right) = {\psi _l \left(\psi_0\right)}, \\
&{\psi _{c,r}\left(\psi_0\right)}\left( {1 + b/2} \right) = {\psi _r \left(\psi_0\right)},
\end{align}
where $\psi_l\left(\psi_0\right)$ and $\psi_r\left(\psi_0\right)$ indicate the left edge and the right edge of the beam with beam focus angle $\psi_0$ when there is no beam squint, i.e., $b=0$. Set $b=0$ in (\ref{beam-range-definition}), and we have
\begin{align}
{\psi _l \left(\psi_0\right)}=\psi_0- \frac{\Delta\psi_{3\text{dB}}}{2},\\
{\psi _r \left(\psi_0\right)}=\psi_0+ \frac{\Delta\psi_{3\text{dB}}}{2},
\end{align}
we have
\begin{align}
&{\psi _{c,l}\left(\psi_0\right)}= \frac{\psi_0- \frac{\Delta\psi_{3\text{dB}}}{2}}{ {1 - \frac{b}{2} } }, \label{psi-cl-calculation} \\
&{\psi _{c,r}\left(\psi_0\right)}= \frac{\psi_0+ \frac{\Delta\psi_{3\text{dB}}}{2}}{{1 +\frac{b}{2}  }  }, \label{psi-cr-calculation} \\
&\psi_0={\left( {1 -\frac{b}{2}} \right)} {\psi _{c,l}} + \frac{\Delta\psi_{3\text{dB}}}{2}.
\label{psi0-calculation}
\end{align}

Similarly, for $\psi_c<0$, we can obtain
\begin{align}
&{\psi _{c,l}\left(\psi_0\right)}= \frac{\psi_0- \frac{\Delta\psi_{3\text{dB}}}{2}}{ {1 +\frac{b}{2} } }, \label{psi-cl-calculation-psic-less-0} \\
&{\psi _{c,r}\left(\psi_0\right)}= \frac{\psi_0+ \frac{\Delta\psi_{3\text{dB}}}{2}}{{1 -\frac{b}{2}  }  }. \label{psi-cr-calculation-psic-less-0}
\end{align}

In considering effect of beam squint, the effective beamwidth for a beam with focus angle $\psi_0$ is
\begin{align}
\Delta\psi_{3\text{dB}} \left(\psi_0\right)&={\psi _{c,r}\left(\psi_0\right)}-{\psi _{c,l}\left(\psi_0\right)} \nonumber \\
&= \left\{
    {\begin{array}{{l}}
    {\frac{{\Delta {\psi _{3{\rm{dB}}}} - b{\left|\psi _0\right|}}}{{1 - \frac{{{b^2}}}{4}}}, \quad {\psi _l}\left( {{\psi _0}} \right){\psi _r}\left( {{\psi _0}} \right) > 0}  \\
    {\frac{{\Delta {\psi _{3{\rm{dB}}}}}}{{1 + \frac{b}{2}}},\quad \quad \quad \ {\psi _l}\left( {{\psi _0}} \right){\psi _r}\left( {{\psi _0}} \right) \le 0,}  \\
    \end{array}}
    \right.
\label{effective-beam-width}
\end{align}
where $\Delta\psi_{3\text{dB}}$ is calculated with (\ref{psi-cl-calculation}) and (\ref{psi-cr-calculation}), or (\ref{psi-cl-calculation-psic-less-0}) and (\ref{psi-cr-calculation-psic-less-0}) if ${\psi _l}\left( {{\psi _0}} \right),{\psi _r}\left( {{\psi _0}} \right) > 0$; otherwise, $\Delta\psi_{3\text{dB}}$ is calculated with (\ref{psi-cr-calculation}) and (\ref{psi-cl-calculation-psic-less-0}). If ${\psi _l}\left( {{\psi _0}} \right),{\psi _r}\left( {{\psi _0}} \right) > 0$, $\Delta\psi_{3\text{dB}} \left(\psi_0\right)$ is not fixed, and it is a linear function of $\psi_0$ for fixed $b$.

Since $\Delta\psi_{3\text{dB}} \left(\psi_0\right) >0$ for all $\psi_0 \in \left[ -\psi_m, \psi_m \right]$, based on (\ref{beam-width-N}) and (\ref{effective-beam-width}), we have the following theorem.
\begin{mytheorem}
If the number of antennas in the array $N$ is fixed, the fractional bandwidth $b$ is upper bounded by
\begin{equation}
b\le \frac{{\Delta\psi_{3\text{dB}}} }{\psi_m }=\frac{1.772}{\psi_m N};
\label{b-upper-bounded}
\end{equation}
Conversely, if the fractional bandwidth $b$ is fixed, the number of antennas in the array $N$ is upper bounded by
\begin{equation}
N \le \floor[\bigg]{\frac{1.772}{\psi_m b}},
\label{N-upper-bounded}
\end{equation}
where $\floor[\big]{\cdot}$ is a floor function.
\label{theorem:b-N-upper-bounded}
\end{mytheorem}

\subsection{Codebook Design with Beam Squint}

\begin{algorithm}[t]
\caption{Codebook Design with Beam Squint}
\label{codebook-design-beam-squint}
\begin{algorithmic}[1]
\Procedure{1}{Odd-Number Codebook}
\State Let $k=1$, $\psi_0=0$
\State Align the first beam so that $\psi_0=0$
\State Calculate $\psi _{c,r}$ based on (\ref{psi-cr-calculation})
\While {$\psi _{c,r}<\psi_m$}
    \State Let $\psi _{c,l}=\psi _{c,r}$
    \State Calculate $\psi_0$ based on (\ref{psi0-calculation})
    \State Calculate $\psi _{c,r}$ based on (\ref{psi-cr-calculation})
    \If {$\psi _{c,l} \ge \psi _{c,r}$}
        \State Codebook does not exist; algorithm ends
    \EndIf
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $\psi_0$
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $-\psi_0$
\EndWhile
\State $\left| \mathcal{C}_{o} \right|=k$
\EndProcedure

\Procedure{2}{Even-Number Codebook}
\State Let $k=0$, $\psi _{c,r}=0$
\While {$\psi _{c,r}<\psi_m$}
    \State Let $\psi _{c,l}=\psi _{c,r}$
    \State Calculate $\psi_0$ based on (\ref{psi0-calculation})
    \State Calculate $\psi _{c,r}$ based on (\ref{psi-cr-calculation})
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $\psi_0$
    \State $k=k+1$
    \State Align the $k$th beam with focus angle $-\psi_0$
\EndWhile
\State $\left| \mathcal{C}_{e} \right|=k$
\EndProcedure
\State ${\left| \mathcal{C} \right| }_{min}=\min\left(\left| \mathcal{C}_{o} \right| , \left| \mathcal{C}_{e} \right| \right)$
\State Select the Procedure that achieves ${\left| \mathcal{C} \right| }_{min}$
\end{algorithmic}
\end{algorithm}

\begin{figure}
\begin{center}
\includegraphics[width=88 mm]{figs/Codebook_design_plot_without_squint}
\end{center}
\caption{Example of a codebook design for the minimum codebook size without consideration of beam squint according to Algorithm~\ref{codebook-design-no-beam-squint}. The two beams with the largest $\psi_0$ are located closely. ${\left| \mathcal{C} \right| }_{min}=19$. $N=16$, $\psi_m=1$, $d=\lambda_c/2$.}
\label{fig:codebook-design-without-wideband-effect}
\bigskip
\begin{center}
\includegraphics[width=88 mm]{figs/Codebook_design_plot_with_squint_even}
\end{center}
\caption{Example of a codebook design for the minimum codebook size with consideration of beam squint according to the Algorithm~\ref{codebook-design-beam-squint}. The beams in the figure are for the carrier frequency. ${\left| \mathcal{C} \right| }_{min}=22$. $N=16$, $\psi_m=1$, $d=\lambda_c/2$. $b=0.0342$ corresponds to $B=$2.5 GHz at $f_c=$73 GHz.}
\label{fig:codebook-design-example}
\end{figure}

Again, to achieve the minimum codebook size, the beam coverages in the codebook should have as little overlap as possible. Algorithm~\ref{codebook-design-beam-squint} shows one method to design a codebook that achieves minimum codebook size with beams symmetric to broadside. Before running the algorithm, it is not clear whether the codebook contains an odd or an even number of beams that achieve the minimum codebook size. Procedure~1 in Algorithm~\ref{codebook-design-beam-squint} designs an odd-number codebook, denoted as ${\mathcal{C}_{o} }$; Procedure~2 designs an even-number codebook, denoted as ${\mathcal{C}_{e} }$. For aligning the beams on the right of the broadside, ${\psi _{c,r}\left(\psi_0\right)}$ of the beam with smaller beam focus angle is equal to ${\psi _{c,l}\left(\psi_0\right)}$ of the beam with larger beam focus angle. The first determined beam in Algorithm~\ref{codebook-design-beam-squint} has $\psi_0=0$.

Either ${\mathcal{C}_{o} }$ or ${\mathcal{C}_{e}}$ achieves the minimum codebook size. The minimum codebook size is also calculated in Algorithm~\ref{codebook-design-beam-squint}, i.e.,
\begin{align}
{\left| \mathcal{C} \right| }_{min}=\min\left(\left| \mathcal{C}_{o} \right| , \left| \mathcal{C}_{e} \right| \right).
\end{align}

\section{Numerical Result}

\begin{figure}[t]
\begin{center}
\includegraphics[width=88 mm]{figs/Codebook_size_calculation_LinearArray_plot_Vs_b_multi_N}
\end{center}
\caption{Minimum codebook size $\left| \mathcal{C} \right|_{min}$ required to cover $\psi \in \left[ -1, 1 \right]$ as a function of fractional bandwidth $b$. The dashed lines indicate upper bounds. $g_t$ is set to be 3~dB smaller than $g_m$. $d=\lambda_c/2$.}
\label{fig:codebook-size-vs-b}
\end{figure}

Fig.~\ref{fig:codebook-design-without-wideband-effect} and \ref{fig:codebook-design-example} show examples of codebook designs for the minimum codebook size without and with considering beam squint, respectively, for $N=16$, $\psi_m=1$ and $d=\lambda_c/2$. $\psi_m=1$ denotes that the maximum targeted covered AoA/AoD of the codebook is $90^{\circ}$. If beam squint is considered, the minimum codebook size increases. The beams are denser for higher focus angles, because of the larger reduction of effective beamwidth for a larger beam focus angle as shown in (\ref{effective-beam-width}). In the two figures, the fractional bandwidth $b$ corresponds to $B=2.5$~GHz and $f_c=73$~GHz. The minimum codebook size in the examples for a $N=16$ ULA is increased by 15.8\% compared to that without considering beam squint, and it is increased by 54\% for a $N=32$ ULA with the same $b$.

Fig.~\ref{fig:codebook-size-vs-b} illustrates examples of the minimum codebook size $\left| \mathcal{C} \right|_{min}$ as a function of fractional bandwidth $b$ for $\psi_m=1$, i.e., the maximum targeted covered AoA/AoD of the codebook is $90^{\circ}$.
$b=0$ is equivalent to not considering beam squint.  For all $N$, as $b$ increases, $\left| \mathcal{C} \right|_{min}$ increases, and $b$ is upper bounded for fixed $N$, confirming Theorem~\ref{theorem:b-N-upper-bounded}. As $b$ approaches to its unreachable upper bound $\frac{\Delta\psi_{3\text{dB}}}{\psi_m}$, $\left| \mathcal{C} \right|_{min}$ increases significantly. If $b \to \frac{\Delta\psi_{3\text{dB}}}{\psi_m}$, the effective beamwidth of the beam with $\psi_0=\psi_m$, $\Delta\psi_{3\text{dB}} \left(\psi_0\right) \to 0$; therefore, $\left| \mathcal{C} \right|_{min}{ \to \infty }$. If $b \ge \frac{\Delta\psi_{3\text{dB}}}{\psi_m}$, no codebook exists.

\begin{figure}
\begin{center}
\includegraphics[width=88 mm]{figs/Codebook_size_calculation_LinearArray_plot_Vs_N_multi_b}
\end{center}
\caption{Minimum codebook size $\left| \mathcal{C} \right|_{min}$ required to cover $\psi \in \left[ -1, 1 \right]$ as a function of the number of antennas $N$. The dashed lines indicate upper bounds. $g_t$ is set to be 3~dB smaller than $g_m$. $d=\lambda_c/2$. $b=0.0179$ corresponds to $B=0.5$ GHz and $f_c=28$ GHz; $b=0.0342$ corresponds to $B=2.5$ GHz and $f_c=73$ GHz; $b=0.0360$ corresponds to $B=2.16$ GHz and $f_c=60$ GHz; $b=0.0714$ corresponds to $B=2$ GHz and $f_c=28$ GHz.}
\label{fig:codebook-size-vs-N}
\end{figure}

Fig.~\ref{fig:codebook-size-vs-N} shows examples of the minimum codebook size $\left| \mathcal{C} \right|_{min}$ as a function of the number of antennas $N$ for $\psi_m=1$, i.e., the maximum targeted covered AoA/AoD of the codebook is $90^{\circ}$. $\left| \mathcal{C} \right|_{min}$ increases as $N$ increases. $N$ is upper bounded by $\floor[\big]{\frac{1.772}{\psi_m b}}$ as described in (\ref{N-upper-bounded}). When $N$ goes beyond the upper bound, no codebook exists. The maximum array gain only depends on $N$. Therefore, the maximum array gain $g_m=\sqrt{N}$ is also upper bounded.

\iftrue
\section{Conclusion}
We analyzed the effect of beam squint in analog beamforming with a phased array implemented via phase shifters as well as its influence on codebook design. Overall, the effect of beam squint increases if the AoA/AoD is away from the beam focus angle, and it tends to increase as AoA/AoD increases. We also developed an algorithm to compensate for the beam squint in codebook design. Beam squint increases the codebook size if the minimum array gain is required for all subcarriers in the wideband signal. The fractional bandwidth is upper bounded for fixed number of antennas; the number of antennas in the phased array is correspondingly upper bounded for fixed fractional bandwidth.

\fi

\balance

\bibliographystyle{IEEEtran}