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

Relativistic nuclear collisions provide a unique opportunity to create and study matter at temperatures and energy densities well above the boundary where hadronic degrees of freedom become irrelevant and are replaced with quark and gluon degrees of freedom. Whereas the energy density inside a hadron might be on the order of 100 MeV/fm$^3$, high-energy nuclear collisions produce mesoscopic regions where the average energy density can exceed 10 GeV/fm$^3$. Unfortunately, they are far from ideal environments for extracting the properties of super-hadronic matter. Only the asymptotic momenta of outgoing particles are experimentally accessible, so three-dimensional models of the dynamics are essential for interpreting results and forming rigorous conclusions about either the nature and properties of the matter or about the evolution itself. One serious limitation is that the exact nature of the initial deposition and dispersion of energy, momentum, charge, and baryon number in heavy-ion collisions is currently only partially understood. Nearly any experimental measurement will be sensitive to these initial conditions as well as to combined contributions from the proceeding stages of rapid expansion and cooling through partially equilibrated QGP and hadronic phases. This entangles contributions from the various components of the underlying physics and obscures the interpretation of many experimental results.

The development of models that can reliably mimic the various complexities of heavy-ion collisions has been a key focus of the field and a necessary step towards condensing experimental results into quantitative statements about the underlying physics. Unfortunately, determining physical quantities like the shear and bulk viscosities of QGP is non-trivial even with high quality models in place. Experimental data can be largely described using these models but the conclusions drawn depend strongly on the validity of the initial state model and other model assumptions. This problem is only confounded by the computational complexity of typical heavy-ion collision models which, in many cases, makes a full exploration of model parameter space prohibitively expensive and, in turn, encourages simplifying assumptions.

The situation of dealing with very complicated models that map high dimensional parameter spaces to rich heterogenous observational data sets is a common one. It's faced in disparate fields of science such as biology, econometrics, and cosmology. The common solution to this problem is to apply a Bayesian inference approach where a posterior distribution over model parameters can be mapped out based on their consistency with a set of observables using a Markov Chain Monte Carlo (MCMC) or a related approach \cite{Habib:2007ca,Gomez:2012ak}. This allows for quantitative conclusions about the underlying parameters to be made given the validity of the model and prior assumptions.

A primary difficulty that has historically limited the applicability of a Bayesian approach in heavy-ion physics has been the large number of model evaluations necessary to accurately map out a high dimensional posterior distribution. This, coupled with the already computationally expensive nature of the simulations, makes a direct application of an MCMC approach essentially impossible with currently available hardware. Recent progress in methodology has overcome this difficulty through the use of Gaussian process emulators which can be trained to accurately reproduce the results of actual simulated models. The reduced computational load has opened the door to quantitative analyses of model parameters which has allowed for systematic constraints of the nuclear equation of state and of the QGP shear viscosity \cite{Novak:2013bqa,Pratt:2015zsa}.

Stating the single point in the multidimensional parameter space that maximizes the likelihood disregards the uncertainty with which the number is stated. In contrast, an MCMC exploration of parameter space provides not only the uncertainties for each parameter, but the covariances with other parameters along with the complete shape of the posterior likelihood. The methods of \cite{Novak:2013bqa} indeed provide the complete likelihood distribution, and are now extracting the field's first rigorous quantitative conclusions from heavy-ion collisions at the highest energies.

Although the Bayesian approach has proved fruitful in nuclear physics, progress in understanding the complexities of generated posterior distributions and how they're driven by the experimental data has been slow. It is relatively straightforward to understand how parameters drive observables through varying them and seeing what changes. Going in the opposite direction, seeing how the measurements drive the inferences about parameters, is more challenging for several reasons. For example, changing a certain parameter might readily alter a given observable, but there may exist some combination of other parameters that can compensate without affecting the remaining observables. Even an observable which does not seem to directly depend on a given parameter, might help constrain other parameters, which then helps constrain the first parameter. In deciding which observables to either measure or to better measure, one would like to know how the width of the posterior likelihood distribution is affected by reducing the uncertainty of a given observable. These relations are critical to gaining insight into understanding not only the degree to which model parameters are being constrained, but how and why.  Here, we present a new set of techniques for addressing all of these questions and then apply the techniques to the 14-parameter analysis first presented in \cite{Pratt:2015zsa}. This analysis represents the field's first quantitative evaluation of much of the field's consensus-based understanding about which observables are truly responsible for addressing the community's most basic questions about the bulk properties of nuclear matter and about the evolution of a high-energy heavy-ion collision.

In the next two sections we review the model, data and techniques used to generate our previous results in \cite{Novak:2013bqa,Pratt:2015zsa}, presenting a wider range of results from the projections of the MCMC procedure than were presented in \cite{Pratt:2015zsa}. In Sec. \ref{sec:sensitivity} we present techniques for determining linear relations between observables and parameters and between uncertainties in observables and widths of the posterior distribution in parameter space. We apply these techniques to the heavy-ion analysis described in Sec.s \ref{sec:model} and \ref{sec:mcmc}, and focus on the ramifications for extracting the equation of state and viscosity. Results are summarized and an outlook is presented in the final section.

\section{Model and Data Overview}
\label{sec:model}

The details of the model are expounded in more detail in \cite{Novak:2013bqa,Pratt:2015zsa}, but we will briefly review the basics and focus on describing the 14 parameters varied in this analysis along with the set of observables. Our model consists of the generation of an initial state which is then fed into a 2-dimensional hydrodynamics simulation followed by transition to a microscopic simulation, known as a hadronic cascade, at a transition temperature of $T_{0}=165\ \text{MeV}$. An assumed symmetry based on an invariance to boosts along the beam axis makes it possible to approximate three-dimensional treatment with a two-dimensional model that discretizes the two transverse coordinates.

Ten model parameters are used to vary the initial state at a time $\tau_0=0.8$ fm/$c$, two are used to determine the equation of state in the hydrodynamic stage, and the last two are used to determine the shear viscosity and it's temperature dependence in the hydrodynamic stage. The shear viscosity is parameterized as
\begin{equation}
\frac{\eta}{s}
=\left(\frac{\eta}{s}\right)_{0}+\eta'\ln\frac{T}{T_0}
\end{equation}
where $(\eta/s)_0$ is the viscosity at $T_0$ and $\eta'$ describes the temperature dependence. These two parameters will be featured in this paper due to their relatively intuitive meaning and relationship with common observables. Two parameters also encapsulate the equation of state. The speed of sound, $c_s$, as a function of the energy density, $\epsilon$, is described as
\begin{eqnarray}
c_s^2(\epsilon)&=&c_s^2(\epsilon_h)+\left(\frac{1}{3}-c_s^2(\epsilon_h)\right)
\frac{X_0x+x^2}{X_0x+x^2+X'^2},\\
\nonumber
X_0&=&X'Rc_s(\epsilon)\sqrt{12},~~x\equiv \ln(\epsilon/\epsilon_h),
\end{eqnarray}
where $\epsilon_h$ and $c_s(\epsilon_h)$ are the energy density and speed of sound of a hadron gas at the transition temperature, $T_0=165$ MeV. These quantities are calculated by considering a gas of non-interacting hadrons using the masses and spins of particles from the Particle Data Group \cite{pdg}. All resonances with masses below 2 GeV/$c^2$ were included. With this prescription, the equation of state is continuous at $T_0$. The parameter $R$ describe the behavior of $c_s$ just above $T_0$ and the parameter $X'$ provides a scale at which the speed of sound approaches $1/3$. Increasing $X'$ lowers the speed of sound mainly at high energy density or temperature, and increasing $R$ increases the speed of sound, mainly just above $T_0$.

Five of the ten initial state parameters apply only to the description of RHIC (Relativistic Heavy Ion Collider at Brookhaven National Laboratory) data. Only data from Au+Au collisions at $100A$ GeV + $100A$ GeV is considered. The other five described the initial state for Pb+Pb collisions at $1.38A$ TeV + $1.38A$ TeV from the LHC (Large Hadron Collider at CERN). The initial transverse energy density profile that instantiates the hydrodynamics has the form \cite{Novak:2013bqa}
\begin{eqnarray}
\epsilon(x,y)&=&f_{\rm wn}\epsilon_{\rm wn}(x,y) + (1-f_{\rm wn})\epsilon_{\rm sat}(x,y),
\end{eqnarray}
where the parameter $f_{\rm wn}$ describes the weighting between two other parameterized forms, $\epsilon_{\rm wn}$ which is the wounded nucleon form \cite{Miller:2007ri}, and $\epsilon_{\rm sat}$, which describes a form more in line with some ideas of saturation, similarly to \cite{Drescher:2007cd}.
\begin{eqnarray}
\label{eq:wn}
\epsilon_{\rm wn}(x,y)&=&
Z_\epsilon\frac{(dE_\perp/dy)_{pp}\sigma_{\rm nn}}{2\sigma_{\rm sat.}}
T_A(x,y)\left(1-\exp(-T_B(x,y)\sigma_{\rm sat})\right),\\
\nonumber
&&+Z_\epsilon\frac{(dE/dy)_{pp}\sigma_{\rm nn}}{2\sigma_{\rm sat.}}T_B(x,y)\left(1-\exp(-T_A(x,y)\sigma_{\rm sat})\right),\\
\label{eq:sat}
\epsilon_{\rm sat}(x,y)&=&Z_\epsilon\frac{(dE/dy)_{pp}\sigma_{\rm nn}}{\sigma_{\rm sat.}}
T_{\rm min}(x,y)\left(1-\exp(-T_{\rm max}(x,y)\sigma_{\rm sat})\right),\\
\nonumber
T_{\rm min}&\equiv&\frac{2T_AT_B}{T_A+T_B},~~T_{\rm max}\equiv(T_A+T_B)/2.
\end{eqnarray}
Here $T_A$ and $T_B$ are the areal densities, the projections of the baryon density of a Au or Pb nucleus onto the transverse plane, and have dimensions of number per area. In the wounded nucleon model, each nucleon that comes within the nucleon-nucleon cross section, $\sigma_{nn}=42$ mb at RHIC energies and 73 mb at LHC energies, of one of the other nucleons, known as a participant, contributes to the energy density. If the parameter $\sigma_{\rm sat}$ is set equal to the nucleon-nucleon cross section then each nucleon can only contribute once and additional collisions do not increase the energy density. Relaxing $\sigma_{\rm sat}<\sigma_{nn}$ allows the particle to contribute multiple times. As $\sigma_{\rm sat}$ approaches zero the form gives binary scaling and the resulting energy density is proportional to $T_AT_B$. In the saturated form, the energy density is principally determined by the smaller of the two areal densities. Thus, if one nucleon overlaps 5 other nucleons, the energy density will be only slightly less than if it overlapped 10.

The parameter $Z_\epsilon$ represents the energy per unit rapidity per nucleon collision in a dilute reaction relative to the measurement of a $pp$ collision, and both forms satisfy the constraint that for a diffuse overlap of nucleons, i.e. $T_A,T_B\ll \sigma_{pp}$ the energy density scales as binary collisions, $\epsilon(x,y)\rightarrow \sigma_{\rm pp}(dE/dy)_{\rm pp}T_AT_B$ with $Z_\epsilon\approx 1$. However, because this is the energy density at $\tau_0=0.8$ fm/$c$, and not the final energy, and because the energy density includes longitudinal motion, it differs from the measured transverse energy,  and $Z_\epsilon$ was allowed to vary over a small range near unity. For our purposes $(dE/dy)_{\rm pp}$ was set to 2.69 and 6.0 GeV for the two beam energies, though this number is somewhat arbitrary because it is multiplied by $Z_\epsilon$. Finally, the fifth parameter describes the initial transverse flow, which is assumed to have the form
\begin{equation}
\label{eq:univflow}
\frac{T_{0i}}{T_{00}}=F_0\frac{- \partial_i T_{00}}{2T_{00}} \tau,
\end{equation}
where $T_{\alpha\beta}$ is the stress-energy tensor. For a traceless stress energy tensor, which would be expected for non-interacting gluon fields or collision-less massless particles or for conformal hydrodynamics, and also assuming boost invariance, $F_0$ would go to unity \cite{Vredevoogd:2008id}. The final parameter adjusted the initial anisotropy of the stress-energy tensor,
\begin{equation}
T_{xx}=T_{yy}=P(1+2\tau'_{xx}), ~~~T_{zz}=P(1-\tau'_{xx}).
\end{equation}
This form is traceless $(1/3)\sum_{i=1,2,3}T_{ii}=P$ and describes the initial shear, which in the Israel Stewart form of hydrodynamics is a dynamical variable which relaxes toward the Navier Stokes value. Parameters are summarized in Table \ref{table:parameters}.

\begin{table}
\begin{tabular}{|c|c|l|}\hline
Parameter & Min, Max (RHIC/LHC values)& Description\\ \hline
$R$ & -0.9, 2.0 & Equation of State \\
$X'$ & 0.5, 5.0 & Equation of State \\
$(\eta/s)_0$ & 0.02, 0.5 & Viscosity at $T_0$\\
$\eta'$ & 0.0, 3.0 & Temperature dependence of viscosity\\
$Z_\epsilon$ & 0.8 / 0.7, 1.25 / 1.4 & Energy normalization ratio\\
$\sigma_{\rm sat}$(mb) & 22 / 38, 44 / 76 & Saturation cross section\\
$f_{\rm wn}$ & 0.0 / 0.0, 1.0 / 1.0 & Weight of wounded-nucleon parameterization\\
$F_0$ & 0.2 / 0.2, 1.0 / 1.0 & Initial transverse flow\\
$\tau'_{xx}$ & 0.0 / 0.0, 1.0 / 1.0 & Initial anisotropy of stress-energy tensor\\ \hline
\end{tabular}
\caption{\label{table:parameters}
Fourteen model parameters}
\end{table}

The model is used to produce a set of 30 observables that correspond to experimental measurements performed at RHIC and the LHC. These encompass distilled information describing spectra, elliptic flow, and femtoscopic correlations in central and mid-central collisions. These observables were chosen because they all have been shown to characterize features of thermalized bulk matter. Fifteen of the observables correspond to Au+Au collisions at RHIC while the other fifteen correspond to Pb+Pb collisions at the LHC. Observables were taken from two centrality classes, the top $0-5\%$ centrality cut and the set of collisions in the $20-30\%$ centrality cut. It was felt that additional data with centrality between these values would be redundant because of the smooth behavior of the observations with centrality. Further,  it would be difficult to assign uncertainties for observables from more peripheral collisions due to the difficulty in justifying hydrodynamic treatments for systems whose size is not larger than a thermal wavelength, and where a significant fraction of the transverse profile is in the corona, which does not fully thermalize.

Three classes of measurements were considered. The first and most basic is spectra and yields. These were distilled to four numbers for each centrality. Three were the mean transverse momenta, $p_t$, for pions, kaons and protons. The average was taken over a finite range of $p_t$, limited at the low end by experimental constraints and cut off at the high end to minimize the effects of jets which are non-thermal features outside of this analysis. The fourth measurement was the multiplicity of pions within the same $p_t$ range. Because chemical equilibrium was assumed, rather than parameterized, and because baryon annihilation was not included, the yields of kaons and protons were not included. In \cite{Novak:2013bqa} it was shown that model runs that had the same $\langle p_t\rangle$ had the same spectral shapes, so no resolving power was lost by considering only one number to characterize the spectra. Those spectra with the same $\langle p_t\rangle$ as the data also provided good descriptions of the experimental spectral shapes \cite{Novak:2013bqa,Pratt:2015zsa}.

The second class of observables is comprised of femtoscopic radii extracted from two-pion correlations at small relative momentum. Extracting these radii has long been a staple of the field \cite{Lisa:2005dd}. The gaussian radii are functions of momentum, and describe the size and shape of the outgoing phase space cloud of the given momentum. The three radii, $R_{\rm out}, R_{\rm side}$ and $R_{\rm long}$, describe the transverse dimension parallel to the momentum, the transverse size perpendicular to the momentum and the size along the beam axis respectively. For this analysis, the radii were averaged over the various bins in transverse momentum, so that three observables encapsulated the data for each centrality class.

The final observable was $v_2$, which characterizes the anisotropic transverse flow which is driven by the elliptical shape of the initial transverse profile. The anisotropy, $v_2\equiv\langle \cos 2\phi\rangle$ in its simplest definition, was only considered in the 20-30\% centrality bin. This bin was chosen because the model used smooth profiles, generated from the average aerial profiles for a given impact parameter, and neglected the lumpy conditions which one would expect from the finite number of scatterers in a Au or Pb nucleus. Even though this bin is probably the least affected by lumps, the experimental value was reduced by a factor of 0.91 to more fairly compare to the model \cite{Pratt:2015zsa}. The next major analysis planned by this group will incorporate fluctuating initial conditions, and in addition to more accurately calculating $v_2$, could also address the higher components $v_3, v_4\cdots$, which are purely driven by the fluctuations. Because $v_2$ rises nearly linearly with transverse momentum, the value of $v_2$ was averaged over $p_t$, with the values weighted by $p_t$, so that bins with higher $p_t$ were relatively more important. The linear weighting was motivated by performing a principal component analysis of the binned values to find the combination that best captured the variability in the model runs. Table \ref{table:pcaobservables} summarizes the 30 observables in this analysis.

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|}\hline
observable & exp. value & $p_t$ weighting & centrality & collaboration\\ \hline
$v_{2,\pi^+\pi^-}$ & 8.14\%& ave. over 11 $p_t$ bins from 160 MeV/$c$ to 1 GeV/$c$ & 20-30\% & STAR \cite{Adams:2004bi}\\
$R_{\rm out}$ & 5.28 fm &ave. over 4 $p_t$ bins from 150-500 MeV/$c$ & 0-5\% & STAR \cite{Abelev:2009tp}\\
$R_{\rm side}$ & 4.81 fm & ave. over 4 $p_t$ bins from 150-500 MeV/$c$& 0-5\% & STAR \cite{Abelev:2009tp} \\
$R_{\rm long}$ & 5.47 fm & ave. over 4 $p_t$ bins from 150-500 MeV/$c$& 0-5\% & STAR \cite{Abelev:2009tp} \\
$R_{\rm out}$ & 4.27 fm & ave. over 4 $p_t$ bins from 150-500 MeV/$c$& 20-30\% & STAR \cite{Abelev:2009tp} \\
$R_{\rm side}$ & 3.99 fm &  ave. over 4 $p_t$ bins from 150-500 MeV/$c$& 20-30\% & STAR \cite{Abelev:2009tp} \\
$R_{\rm long}$ & 4.53 fm & ave. over 4 $p_t$ bins from 150-500 MeV/$c$& 20-30\% & STAR \cite{Abelev:2009tp} \\
$\langle p_t\rangle_{\pi^+\pi^-}$ & 494.4 MeV & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 0-5\% & PHENIX \cite{Adler:2003cb} \\
$\langle p_t\rangle_{K^+K^-}$ & 796 MeV & 0.4 GeV/$c < p_t< 1.6$ GeV/$c$ & 0-5\% & PHENIX \cite{Adler:2003cb} \\
$\langle p_t\rangle_{p\bar{p}}$ & 1.135 GeV & 0.6 GeV/$c < p_t< 2.0$ GeV/$c$ & 0-5\% & PHENIX \cite{Adler:2003cb} \\
$\langle p_t\rangle_{\pi^+\pi^-}$ & 487.5 MeV & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 20-30\% & PHENIX \cite{Adler:2003cb} \\
$\langle p_t\rangle_{K^+K^-}$ & 792 MeV &0.4 GeV/$c < p_t< 1.6$ GeV/$c$ & 20-30\% & PHENIX \cite{Adler:2003cb} \\
$\langle p_t\rangle_{p\bar{p}}$ & 1.111 GeV & 0.6 GeV/$c < p_t< 2.0$ GeV/$c$ & 20-30\% & PHENIX \cite{Adler:2003cb} \\
$\pi^+\pi^-$ yield & 422 & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 0-5\% & PHENIX \cite{Adler:2003cb}\\
$\pi^+\pi^-$ yield & 188.7 & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 20-30\% & PHENIX \cite{Adler:2003cb}\\
\hline
$v_{2,\pi^+\pi^-}$ & 9.56\% & ave. over 11 $p_t$ bins from 0.15 to 1.2 GeV/$c$ & 20-30\% & ALICE \cite{Abelev:2014pua}\\
$R_{\rm out}$ & 5.46 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$ & 0-5\% & ALICE \cite{Graczykowski:2014hoa} \\
$R_{\rm side}$ & 5.32 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$  & 0-5\% & ALICE \cite{Graczykowski:2014hoa}  \\
$R_{\rm long}$ & 5.72 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$  & 0-5\% & ALICE \cite{Graczykowski:2014hoa}  \\
$R_{\rm out}$ & 4.07 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$  & 20-30\% & ALICE \cite{Graczykowski:2014hoa}  \\
$R_{\rm side}$ & 4.12 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$  & 20-30\% & ALICE \cite{Graczykowski:2014hoa}  \\
$R_{\rm long}$ & 4.41 fm & ave. over 7 $p_t$ bins from 200 to 900 MeV/$c$  & 20-30\% & ALICE \cite{Abelev:2009tp} \\
$\langle p_t\rangle_{\pi^+\pi^-}$ & 459.1 MeV & 0.1 GeV/$c < p_t< 1.2$ GeV/$c$ & 0-5\% & ALICE \cite{Abelev:2013vea} \\
$\langle p_t\rangle_{K^+K^-}$ & 775 MeV & 0.2 GeV/$c < p_t< 1.6$ GeV/$c$ & 0-5\% & ALICE \cite{Abelev:2013vea} \\
$\langle p_t\rangle_{p\bar{p}}$ & 1.137 GeV & 0.2 GeV/$c < p_t< 2.0$ GeV/$c$ & 0-5\% & ALICE \cite{Abelev:2013vea} \\
$\langle p_t\rangle_{\pi^+\pi^-}$ & 455 MeV & 0.1 GeV/$c < p_t< 1.2$ GeV/$c$ & 20-30\% & ALICE \cite{Abelev:2013vea} \\
$\langle p_t\rangle_{K^+K^-}$ & 758 MeV & 0.2 GeV/$c < p_t< 1.6$ GeV/$c$ & 20-30\% & ALICE \cite{Abelev:2013vea} \\
$\langle p_t\rangle_{p\bar{p}}$ & 1.110 MeV & 0.2 GeV/$c < p_t< 2.0$ GeV/$c$ & 20-30\% & ALICE \cite{Abelev:2013vea} \\
$\pi^+\pi^-$ yield & 1258 & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 0-5\% & ALICE \cite{Abelev:2013vea}\\
$\pi^+\pi^-$ yield & 523 & 0.2 GeV/$c < p_t< 1.2$ GeV/$c$ & 20-30\% & ALICE \cite{Abelev:2013vea}\\ \hline
\end{tabular}
\caption{\label{table:pcaobservables} Observables used to compare models to data from RHIC (STAR and PHENIX collaborations) and from ALICE (at the LHC). To account for non-flow correlations, the value of $v_2$ was reduced by 9\% from the experimental values to account for the approximation of smooth initial conditions.}
\end{table}

\section{Markov Chain Monte Carlo Procedure and Results}
\label{sec:mcmc}

The standard approach to determining the likely regions of parameters, $\vec{x}$, from comparing model values, $\vec{y}_M(\vec{x})$, to experimental values, $\vec{y}_{\rm exp}$, is through MCMC, which provides a sampling of parameters $\vec{x}$ that are chosen weighted proportional to the likelihood. For this study the likelihood is chosen to have a simple Gaussian form,
\begin{equation}
{\mathcal L}(\vec{x})\sim \exp\left\{-\sum_a\frac{(y_{M,a}(\vec{x})-y_{{\rm exp},a})^2}{2\sigma_a^2}\right\}.
\end{equation}
The uncertainties incorporate both experimental uncertainties and the shortcomings of the model. For instance, if the equation of state and viscosity were perfectly described, and if the initial state was parameterized most correctly, the model uncertainty describes how accurately one would expect to reproduce the final-state observables given the missing physics in the model and the uncertainties in the measurements. In this study, which uses the same model output as in \cite{Pratt:2015zsa}, the uncertainties $\sigma_a$ were all set to 6\% of the experimental value. A more accurate determination of the uncertainty would require detailed model studies to estimate the impact of missing physics, and a more detailed understanding of experimental uncertainties. It would not be surprising to find that a few of these observables are understood with slightly better uncertainty or that some are somewhat more uncertain. Such an improvement should involve input from both the experimental and modeling communities. In \cite{Novak:2013bqa} the analysis was repeated with 9\% uncertainties and the widths of the resulting posterior distributions only increased by $\approx 20\%$.

Due to the infeasibility of performing millions of full model runs, the MCMC procedure implemented an emulator in the place of the full model \cite{Novak:2013bqa,Gomez:2012ak}. The emulator uses an interpolation algorithm to determine the observables from 1200 full-model runs. The first 1000 parameter sets for the full-model runs uniformly covered the 14-dimensional space according to latin hyper-cube sampling. The last 200 parameter sets were chosen to be consistent with the likelihood as calculated from the earlier runs and to provide better coverage in the likely region. To improve efficiency, the emulator calculates principal components rather than each observable. In this way, those linear combinations of observables that stay constant throughout the model runs can be neglected. Thus, rather than evaluating all 30 observables, only 14 principal components were analyzed for this study. Even that number could be reduced because the last several components only varied by a few percent of one value of the experimental uncertainty.

\begin{figure}
\centerline{\includegraphics[width=\textwidth]{ll}}
\caption{\label{fig:ll}
One- and two-dimensional posterior likelihood projections from the MCMC procedure. The one-dimensional projections (along the diagonal) illustrate the degree to which an individual parameter is constrained by comparing to data, whereas the off-diagonal elements illustrate how some linear combinations of two parameters might be either poorly or well constrained. The colored boundaries delineate one-, two- and three-sigma regions, where $n$-sigma refers to likelihoods of $e^{-n^2/2}$ of the maximum likelihood.}
\end{figure}
The MCMC procedure was performed with several million random steps according to the Metropolis algorithm. This provides a sampling of the posterior distribution displayed in Fig. \ref{fig:ll}. In addition to projections onto one dimension of the parameter space, the off-diagonal plots show two-dimensional projections. When the elliptic shapes of the two-dimensional projections lie at an angle, it shows that a specific linear combination of parameters may be well constrained, whereas the orthogonal combination may be poorly constrained. For example, the region of high likelihood for the projection onto the plane of two equation of state parameters, one can see that if one increases $R$ and $X'$ by similar percentages that the likelihood changes little.  Similarly, for the two viscosity parameters one finds that one can find a good match with $(\eta/s)_0\approx 0.2$ with a modest temperature dependence, $\eta'$. One can also match with lower values of $(\eta/s)_0$ if the viscosity then rises with temperature. Another instance where the posterior is off-diagonal is in the projection of $(\eta/s)_0$ and $f_{\rm wn}$. The resulting covariance corroborates the arguments that were put forward in \cite{Drescher:2007cd}.

\section{Sensitivity Studies}
\label{sec:sensitivity}

The principal goal of this paper is to present various methods for understanding the role certain observables, or sets of observables, play in constraining the posterior likelihood of given parameters, or sets of parameters. The most straight-forward method to determine the sensitivity is to perform the analysis both with and without a given observable, or set of observables. An example of this is illustrated in Fig. \ref{fig:ll_RHIC_LHC}. The two-dimensional posterior likelihood projection for the two viscosity parameters, $(\eta/s)_0$ and $\eta'$, are redisplayed along with projections where either RHIC data or LHC data are ignored. This makes it clear that the LHC data is especially important for constraining the temperature dependence of the viscosity, $\eta'$.
\begin{figure}
\centerline{
\includegraphics[width=0.16667\textwidth]{ll_visc_nolhc}
\includegraphics[width=0.16667\textwidth]{ll_visc_norhic}
\includegraphics[width=0.16667\textwidth]{ll_visc}
\includegraphics[width=0.166667\textwidth]{ll_eos_nolhc}
\includegraphics[width=0.166667\textwidth]{ll_eos_norhic}
\includegraphics[width=0.166667\textwidth]{ll_eos}}
\caption{\label{fig:ll_RHIC_LHC}
Posterior likelihood distributions for the shear viscosity parameters generated using RHIC and LHC data separately (a \& b) and together (c). We can see that the viscosity at $T_0$ is well constrained by RHIC data alone but that there is very little constraint of its temperature dependence without the LHC data. This is consistent with the fact that LHC data probe higher temperatures. Combining the datasets clearly constrains the data better than a single data set. Panels (d-e) show corresponding information for the two equation of state parameters. In this case it is clear that LHC data provide the bulk of the resolving power. Again, this is expected because the equation of state is fixed at $T_0=165$ MeV, and LHC collisions probe regimes further into the region where the equation of state differs across the parameter space. One can also see that the high degree of covariance in both figures (c) and (f), especially for the covariance between the two equation-of-state parameters in (f). This shows that even though neither parameter is particularly well constrained, there is a linear combination that is well constrained, while the orthogonal linear combination is poorly constrained.}
\end{figure}

However, the method of repeating the statistical analysis without a given observable can be onerous. If one wishes to study the sensitivity to $N_Y=30$ observables, one would repeat the MCMC procedure $N_Y$ times. Here, we present methods for finding three simple measures of the sensitivity using the output from a single MCMC trace. That information would be of the form of a list of values, $\alpha=1\cdots N_{\rm MCMC}$, where $N_{\rm MCMC}$ might be of the order of one million. For each point one would have the parametes $x_{\alpha,i}$ and the observables $y_{\alpha,a}$.

The first measure would simply describe how a given observable, $y_a$, would change due to a change in a parameter $x_i$, while keeping all other parameters constant. Because $y_{M,a}^{\rm (\vec{x})}$ is not purely linear, one needs to choose over what region the derivative is evaluated. Two choices of interest might be the prior or the posterior. Here, we use single brackets, $\langle\cdots\rangle$, to denote an average of the prior, and double brackets, ${\langle\langle}\cdots{\rangle\rangle}$, to denote averages taken over the posterior. For averages over the prior, one can use the lists from the points at which full-model runs were performed for the purpose of training the emulator. If calculations were performed randomly throughout the prior, one could consider the covariance $\langle\delta y_a\delta x_i\rangle$, where $\delta y_a=y_a-\langle y_a\rangle$ and $\delta x_i=x_i-\langle x_i\rangle$. One can calculate the partial derivative of $y_a$ with respect to any parameter by a simple matrix inversion,
\begin{eqnarray}
\label{eq:dydx}
\langle \delta y_a\delta x_i\rangle&=&\sum_j\left\langle\frac{\partial y_a}{\partial x_j}\right\rangle\langle\delta x_j\delta x_i\rangle\\
\nonumber
\left\langle\frac{\partial y_a}{\partial x_j}\right\rangle&=&\langle \delta y_a\delta x_i\rangle \langle\delta x\delta x\rangle^{-1}_{ij}.
\end{eqnarray}
Here, the brackets around $\partial y/\partial x$ note that this is the best slope for this region of parameters space, which in this case is the prior. This expression matches the usual least-squares expression for the slope of a line in multidimensions. Equation \ref{eq:dydx} can be easily altered to address how the slope looks when focused on the posterior region by making the change $\langle\cdots\rangle\rightarrow\langle\langle\cdots\rangle\rangle$.

\subsection{Response of extracted parameter values to changes in experimental measurement}

Equation \ref{eq:dydx} shows how a specific model value, $y_a$, changes when a specific parameter, $x_i$, is varied. However, describing how a small change in an experimental measurement, $y^{\rm(exp)}_a$, affects the average value of a parameter, $\langle\langle x_i\rangle\rangle$, in the average posterior value is a different question. This involves understanding $\partial\langle\langle x_i\rangle\rangle/\partial y^{\rm(exp)}_a$. Here we consider a general expression for any posterior-averaged function of  $\vec{x}$, $\langle\langle f(\vec{x})\rangle\rangle$, which for this specific consideration will be $f(\vec{x})=x_i$.
\begin{eqnarray}
\label{eq:dxdy_general}
\frac{\partial}{\partial y^{\rm(exp)}_a}
\langle\langle f(\vec{x})\rangle\rangle
&=&\frac{\partial}{\partial y^{\rm(exp)}_a}\frac{\int d\vec{x}~f(\vec{x}){\mathcal L}(\vec{x})}{\int d\vec{x}~{\mathcal L}(\vec{x})}\\
\nonumber
&=&
\langle\langle f(\vec{x})\frac{1}{\mathcal L}\frac{\partial}{\partial y^{\rm(exp)}_a}{\mathcal L}\rangle\rangle
-\langle\langle f(\vec{x})\rangle\rangle\langle\langle \frac{1}{\mathcal L}\frac{\partial}{\partial y^{\rm(exp)}_a}{\mathcal L}\rangle\rangle.
\end{eqnarray}
Both terms on the right-hand-side of Eq. (\ref{eq:dxdy_general}) can be calculated from the trace. By setting $f=\delta x_i=(x_i-\bar{x}_i)$, with $\bar{x}_i\equiv\langle\langle x_i\rangle\rangle$, the second term vanishes in this case, and if one has an expression for the likelihood, one need only average $\delta x_i(1/{\mathcal L})\partial{\mathcal L}$ over the sampling of points from the MCMC trace to calculate the required result. One can invoke the emulator to calculate $y(\vec{x})$, and thereby determine both ${\mathcal L}$ and $\partial{\mathcal L}$.

For Gaussian likelihoods Eq. (\ref{eq:dxdy_general}) becomes especially simple. In that case
\begin{eqnarray}
\nonumber
{\mathcal L}&=&\exp\left\{\frac{1}{2}(y_a-y_a^{\rm(exp)})\Sigma_{ab}^{-1}(y_b-y_b^{\rm(exp)})\right\}\\
\nonumber
\frac{1}{\mathcal L}\frac{\partial}{\partial y^{\rm(exp)}_a}{\mathcal L}&=&
\Sigma^{-1}_{ab}(y_a-y^{\rm(exp)}_a),\\
\label{eq:dxdy_offdiag}
\frac{\partial}{\partial y^{\rm(exp)}_a}\bar{x}_i
&=&\Sigma^{-1}_{ab}\langle\langle \delta x_i\delta y_b\rangle\rangle,\\
\nonumber
&=&\Sigma^{-1}_{ab}\left\langle\left\langle\frac{\partial y_a}{\partial x_j}\right\rangle\right\rangle A_{ji},\\
\nonumber
A_{ij}&\equiv&\langle\langle\delta x_i\delta x_j\rangle\rangle.
\end{eqnarray}
Here, $\delta y$ can be relative to any point because $\langle\langle\delta x\rangle\rangle=0$, and the derivative in the last line is the average slope for the posterior region. Eq. (\ref{eq:dxdy_offdiag}) is straightforward to calculate from the output of the MCMC. For the case where $\Sigma$ is diagonal,
\begin{equation}
\label{eq:dxdy}
\left.
\frac{\partial \bar{x}_i}{\partial y_a^{\rm(exp)}}
\right|_{y^{\rm(exp)}_{b\ne a}}=
\frac{1}{\sigma_a^2}\langle\langle\delta y_a\delta x_i\rangle\rangle
\end{equation}

To better quantify how a specific observable constrains a given parameter, one may multiple the expression for $\partial \bar{x}/\partial y^{\rm(exp)}$ in Eq. (\ref{eq:dxdy}) by the measure of how much $y_a$ changes throughout the prior, $\langle\delta y_a^2\rangle^{1/2}$, because model values of an observable must change within the range of the prior if that observable is to provide resolving power. Figure \ref{fig:dydxdxdy} displays this quantity for all observables. One can also see whether the change is positive or negative, which is not obvious and may have an opposite sign compared to $\partial y/\partial x$.
\begin{figure}
\centerline{\includegraphics[width=\textwidth]{dydxdxdy}}
\caption{\label{fig:dydxdxdy}
Model responses of an observable with respect to a given parameter, $\langle\langle\partial y_a/\partial x_i\rangle\rangle$ are displayed in the left-side panel. In the right-side panel, the change of the inferred value of a parameter with respect to changes in a measurement, $\partial\langle\langle x_i\rangle\rangle/\partial y^{\rm(exp)}_a$, are scaled by the spread of model values throughout the prior, $\langle\delta y_a^2\rangle^{1/2}$. Larger absolute values point to measurements which play important roles in constraining that parameter. The signs of the response are not always equal for the corresponding derivative in the two plots.
}
\end{figure}

As examples of the sensitivity analysis one can look at the elliptic flow measured at RHIC and at the LHC. In Fig. \ref{fig:dydxdxdy} one can see that changing the measurement of $v_2$ at RHIC strongly affects the extracted value of the viscosity at $T_0$, $(\eta/s)_0$, but has little effect on the extracted temperature dependence, $\eta'$. In contrast, the measurement of $v_2$ at the LHC more strongly affects $\eta'$, while playing a more minor role in determining $(\eta/s)_0$. Given the fact that the LHC explores higher energy densities, this was expected. One can also see that the constraints on the viscosity were driven by measurements of $v_2$ and the multiplicities, whereas constraints on the two equation-of-state parameters were driven by a wider variety of measurements.

\subsection{Relation between experimental uncertainties and widths of posterior parameter distributions}

We now consider how the uncertainty of a specified observable affects the width of the posterior parameter distribution, i.e., calculate the derivative
\begin{eqnarray}
\label{eq:Rijab}
\frac{\partial}{\partial\Sigma_{ab}}A_{ij}
&=&\langle\langle \delta x_i\delta x_j \frac{1}{\mathcal L}\frac{\partial{\mathcal L}}{\partial\Sigma_{ab}}\rangle\rangle
-\langle\langle \delta x_i\delta x_j\rangle\rangle\langle\langle\frac{1}{\mathcal L}\frac{\partial{\mathcal L}}{\partial\Sigma_{ab}}\rangle\rangle~,
\end{eqnarray}
where the last step followed the steps from the previous subsection used to derive Eq. (\ref{eq:dxdy_general}). Equation (\ref{eq:Rijab}) is fairly straight-forward to calculate, but can be a bit cumbersome. To find a simpler, though approximate, expression we assume that the prior distribution of parameters is Gaussian  even though the results shown here all assumed hard cutoffs, and that $y(\vec{x})$ is linear, so that the posterior distribution is Gaussian. We also use Gaussian likelihoods for comparing to experimental values. When comparing calculations with and without these approximations, results changed little.

With the approximations, the overall likelihood is then of a Gaussian form,
\begin{eqnarray}
\label{eq:gaussform}
{\mathcal L}\sim \exp\left\{-|\delta x|^2/2R^2-\left\langle\left\langle\frac{\partial y_a}{\partial x_i}\right\rangle\right\rangle\Sigma_{ab}^{-1}
\left\langle\left\langle\frac{\partial y_b}{\partial x_j}\right\rangle\right\rangle
\frac{\delta x_i\delta x_j}{2}\right\},\\
\nonumber
\delta x_i=x_i-\bar{x}_i.
\end{eqnarray}
We have assumed that the parameters $x_i$ have all been scaled to have the same prior width $R$. Further, in the neighborhood of the maximum likelihood we assume that $y$ behaves linearly. The width of the posterior, $A_{ij}$, can be taken from the Gaussian form,
\begin{eqnarray}
\nonumber
A_{ij}&=&
\frac{\int d^Nx~\exp\{-A^{-1}_{ij}\delta x_i\delta x_j/2\}\delta x_i\delta x_j}
{\int d^Nx~
\exp\{-A^{-1}_{ij}\delta x_i\delta x_j/2\}},\\
\label{eq:Adef}
A_{ij}^{-1}&\equiv&\delta_{ij}/R^2+\left\langle\left\langle\frac{\partial y_a}{\partial x_i}\right\rangle\right\rangle\Sigma^{-1}_{ab}
\left\langle\left\langle\frac{\partial y_b}{\partial x_j}\right\rangle\right\rangle.
\end{eqnarray}
Our stated goal is to find an expression describing how $A_{ij}$ responds to changes in $\Sigma_{ab}$,
\begin{eqnarray}
\frac{\partial A}{\partial\Sigma_{ab}^{-1}}&=&\frac{\partial}{\partial\Sigma^{-1}_{ab}}(AA^{-1}A)\\
\nonumber
&=&2\frac{\partial}{\partial\Sigma^{-1}_{ab}}A+A\left(\frac{\partial}{\partial\Sigma_{ab}^{-1}}A^{-1}\right)A\\
\nonumber
&=&-A\left(\frac{\partial}{\partial\Sigma^{-1}_{ab}}A^{-1}\right)A\\
\label{eq:ASigmaresult}
\frac{\partial}{\partial\Sigma^{-1}_{ab}}A_{i\ell}&=&-A_{ij}\left\langle\left\langle\frac{\partial y_a}{\partial x_j}\right\rangle\right\rangle\left\langle\left\langle
\frac{\partial y_b}{\partial x_k}\right\rangle\right\rangle A_{k\ell}.
\end{eqnarray}
One can now insert the expression for $\partial y/\partial x$ given in Eq. (\ref{eq:dydx}) into Eq.s \ref{eq:Adef} and \ref{eq:ASigmaresult} to find the needed result.  When working with the hard cutoffs for priors rather than Gaussians, we scale all the parameters to have the same prior width, $\langle x^2\rangle=R^2$, in addition to being centered at zero, $\langle x\rangle=0$. If one uses the prior distribution to calculate $\partial y/\partial x$, then the response in Eq. (\ref{eq:ASigmaresult}) needn't ever access the experimental value, and one doesn't need to perform the MCMC.

\begin{figure}
\centerline{\includegraphics[width=0.6\textwidth]{sensitivity}}
\caption{\label{fig:sensitivity}
The resolving power of a particular observable in determining a specific parameter, $R_{i;a}$ defined in Eq. (\ref{eq:Riadef}), is displayed for all observables and parameters.
}
\end{figure}
If one wishes to use the posterior distribution to calculate the derivatives, the form for the response simplifies
\begin{eqnarray}
\label{eq:ASigmaresult_posteriorA}
\frac{\partial}{\partial\Sigma^{-1}_{ab}}A_{ij}&=&-\langle\langle \delta x_i\delta y_a\rangle\rangle
\langle\langle\delta y_b \delta x_j\rangle\rangle,
\end{eqnarray}
which can also be transformed into an expression in terms of derivatives with respect to $\Sigma$,
\begin{eqnarray}
\nonumber
\frac{\partial}{\partial\Sigma_{ab}}A_{ij}&=&
\frac{\partial A_{ij}}{\partial\Sigma^{-1}_{cd}}\frac{\partial\Sigma^{-1}_{cd}}{\partial\Sigma_{ab}}\\
\nonumber
\frac{\partial\Sigma^{-1}_{cd}}{\partial\Sigma_{ab}}&=&\frac{\partial}{\partial\Sigma_{ab}}\left(\Sigma^{-1}_{ce}\Sigma_{ef}\Sigma^{-1}_{fd}\right)\\
\nonumber
&=&-\Sigma^{-1}_{ca}\Sigma^{-1}_{bd},\\
\label{eq:ASigmaresult_posterior}
\frac{\partial}{\partial\Sigma_{ab}}A_{ij}&=&\Sigma^{-1}_{ac}\langle\langle \delta y_c\delta x_i\rangle\rangle
\langle\langle\delta y_d \delta x_j\rangle\rangle\Sigma^{-1}_{db}.
\end{eqnarray}
In this analysis the uncertainty matrix is diagonal, $\Sigma_{ab}=\sigma_a^2\delta_{ab}$, and Eq. (\ref{eq:ASigmaresult_posterior}) can be used to calculate the resolving power for determining a single parameter $x_i$ due to a single parameter $y_a$.
\begin{eqnarray}
\label{eq:Riadef}
R_{i;a}&\equiv&\frac{\sigma_a}{\langle(\delta x_i)^2\rangle}\frac{\partial}{\partial\sigma_a}A_{ii},\\
\nonumber
&=&\frac{\sigma_a}{\langle(\delta x_i)^2\rangle}
\Sigma^{-1}_{ac}\langle\langle \delta y_c\delta x_i\rangle\rangle
\langle\langle\delta y_d \delta x_j\rangle\rangle\Sigma^{-1}_{da}.
\end{eqnarray}
Results for $R_{i;a}$ are shown in Fig. \ref{fig:sensitivity}.

Even when a specific parameter is not well constrained, a combination of that parameter with another might still be well constrained.  The two equation-of-state parameters are a good example, and the two viscosity parameters also have a strong off-diagonal component to their likelihood projections. These two instances are shown in panels (c) and (f) in Fig. \ref{fig:ll_RHIC_LHC}. If one then wants to understand the contribution of a given observable in constraining the two-dimensional projection, one needs to define a quantity which effectively measures the likely area of the two-dimensional projection onto the parameters $x_i$ and $x_j$. Here we use the determinant of the two-by-two matrix constructed from the $ij$ components of $A_{ij}$.
\begin{eqnarray}
\label{eq:Dijdef}
D_{ij}&\equiv&A_{ii}A_{jj}-A_{ij}A_{ji},\\
\label{eq:Rijadef}
R_{ij;a}&\equiv&\sigma_a\frac{\partial}{\partial\sigma_a}D_{ij}
=\sigma_a\left\{A_{ii}\frac{\partial}{\partial\sigma_a}A_{jj}
+A_{jj}\frac{\partial}{\partial\sigma_a}A_{ii}
-2A_{ij}\frac{\partial}{\partial\sigma_a}A_{ji}
\right\}.
\end{eqnarray}
$D_{ij}$ would represent the product of the eigenvalues of the two-by-two matrix, or equivalently, the square of the area covered by the projection. The way in which $D_{ij}$ changes with respect to a given uncertainty can then be calculated with the help of Eq. (\ref{eq:ASigmaresult_posterior}). These sensitivities are illustrated in Fig. \ref{fig:eos2eta2}. This clearly demonstrates that the extracted viscosity is strongly affected by the $v_2$ measurements, and that multiplicities are also important. Other observables are of secondary, but not negligible importance. For constraining the equation of state, femtoscopic radii seem to provide the most resolving power, but all other observables contribute significantly. This underscores the importance of simultaneously considering multiple classes of observables to constrain the parameter space.

\begin{figure}
\centerline{\includegraphics[width=0.6\textwidth]{eos2eta2}}
\caption{\label{fig:eos2eta2}
The resolving power for determining a pair of parameters, $x_i,x_j$, due to a specific observable $y_a$ as defined by $R_{ij;a}$ in Eq. (\ref{eq:Rijadef}). The sensitivity for the two viscosity parameters is shown in the upper panel and that for the two equation-of-state parameters is displayed in the lower panels. Information about eh viscosity is most strongly determined by measurements of $v_2$ and the multiplicities. The equation of state is most strongly constrained by femtoscopic radii.}
\end{figure}

\section{Summary and Outlook}
\label{sec:summary}

The results here are unprecedented for this field, and for the first time illustrate a systematic method for identifying the critical links between parameters and observables. Remarkably, the links identified by the procedure reinforced the general knowledge of the field. For instance, it was indeed found that measurements of the elliptic anisotropy provide strong constraints of the viscosity. Further, the expectation that RHIC data would play a larger relative role in determining the viscosity near $T_c$, and that the LHC data would play a larger role in determining the temperature dependence was confirmed. The relatively large role of the multiplicities was not necessarily expected, nor was the fact that other observables provide non-negligible resolving power.

Whereas there was consensus within the field that $v_2$ would be the most important observable to determine the viscosity, the role of various observables in constraining the equation of state was very much in dispute. Figures \ref{fig:sensitivity}, \ref{fig:dydxdxdy} and \ref{fig:eos2eta2} all show that the femtoscopic source sizes provide the most resolving power. This validates ideas based on low pressures leading to more elongated phase space distributions in the outward direction and that the longitudinal sizes would increase if transverse expansion was slower, \cite{Pratt:1986cc,Pratt:2008qv,Lisa:2005dd}. Additionally, source sizes played a role in measuring the final entropy, which is also a measure of the equation of state \cite{Pal:2003rz}. The mean transverse momentum had long ago been pointed out as being sensitive to the temperature, and therefore the equations of state, \cite{VanHove:1982vk}, and even $v_2$ had been suggested as a means for extracting early pressure \cite{Sorge:1996pc}. Therefore, it was not surprising to see the resolving power for the equation of state in Fig. \ref{fig:eos2eta2} spread across all observables.

The techniques presented here could play a pivotal role in determining the direction of future experiments. Before embarking on an expensive experimental program to improve the measure of a specific observable, one could check to see how that improvement might indeed better determine parameters of greatest interest. For example, one could determine whether running an accelerator for a specific projectile target combination, or at a new beam energy, or with higher statistics would best provide insight into better determining the equation of state. In many cases it would behoove the community to pre-analyze a project with modern statistical techniques before investing the cost and manpower for the effort.

\begin{acknowledgments}
This work was supported by the National Science Foundation's Cyber-Enabled Discovery and Innovation Program through grant NSF-0941373 and by the Department of Energy Office of Science through grant number DE-FG02-03ER41259.
\end{acknowledgments}

\input{biblio.tex}

\end{document}


\section{Introduction}
A topological insulating (TI) phase is a new state of quantum matter featuring massless Dirac-like boundary states whose robustness is guaranteed by the time-reversal symmetry (TRS)\cite{Hasan10,Qi11,Moore10,Ando13}. In this rapidly growing field, search for new topological phases has recently turned to extending the consideration of symmetry-protected states to include non-spatial symmetries\cite{Schnyder08,Kitaev09} such as, the particle-hole and chiral symmetries as well as the spatial symmetries\cite{Hughes10,Turner12,Fang12,Slager13,Sato14}. In particular, the nontrivial band topology can be shown to be protected by certain crystal symmetries, leading to the new class of TIs called topological crystalline insulators (TCIs)\cite{LFTCI,mChen,Chiu13}. Like the TIs, a three dimensional (3D) TCI is also predicted to host metallic surface states on surfaces which preserve appropriate crystal symmetries\cite{Hsieh12}.

To date, the semiconducting Pb$_{1-x}$Sn$_x$(Te, Se) alloys are the only experimentally realized TCI materials, whose topological nature has been verified through direct observation of Dirac-like surface states via angle-resolved photoemission spectroscopy (ARPES) experiments\cite{Xu12a,Dziawa12,Tanaka12}. These gapless surface states are protected by the mirror symmetry with respect to the (110) or ($\bar{1}$10) lattice plane. The characteristic features of the associated Dirac states such as linear dispersion, Lifshitz transition, spin/orbital texture, among others, have been examined via scanning tunneling microscopy/spectroscopy (STM/STS) and transport measurements\cite{Xu12a,Yazdani13,Taskin14,OkadaTCIstrain,ZeljkovicTCIphase,ZeljkovicTCIorbital}. Moreover, it has been found that, below a critical temperature, which depends on chemical composition, the cubic lattice structure can be distorted by strains, resulting in an orthorhombic or a rhombohedral structure\cite{Iizumi75}. Remarkably, in a recent STM experiment by Okada {\it et al.}\cite{OkadaTCIstrain}, two unexpected Landau levels have been observed for the (001) surface of Pb$_{1-x}$Sn$_x$Se under a perpendicular magnetic field. The presence of these two extra levels is believed to be associated with the gap opening of the two surface Dirac cones induced by a ferroelectric-like lattice distortion, which breaks the corresponding mirror symmetry. It is clear that it is important to understand effects of strains on the electronic structure of TCIs in order to gain a deeper handle on the nature and origin of their topological states, and how these states could be manipulated for practical applications of the TCIs.

Despite some strain-related studies in the literature\cite{Fang14,Serbyn14}, we are not aware of a systematic investigation of how spin-textures and other properties of Dirac states in the TCIs evolve under various symmetry breaking strains. So motivated, here we examine strain effects on the (001) surface of the Pb$_{1-x}$Sn$_x$(Te, Se) TCIs. Our analysis is based on an effective four-band model, which is shown to capture all essential features of the topological surface states. The form of strain related perturbations on the Hamiltonian is clarified through general symmetry considerations, allowing us to delineate how characteristic features of the topological states in TCIs evolve under strains.

Significantly, our analysis shows that certain mirror-symmetry-breaking strains induce hedgehog spin texture with out-of-the-plane spin-polarization at the Dirac point. Such a spin texture of Dirac states has been reported previously in a 3D TI, but only in the presence of an exchange field\cite{Xu12b}. Our findings thus identify a possible new pathway for realizing TCI-based spintronics devices without requiring a TRS breaking field\cite{Liu14}. Moreover, we show that the interplay of strains and applied Zeeman field can be used to tune the Chern number of the surface ground state, which is a topological invariant characterizing a quantum anomalous Hall (QAH) or a quantum Hall insulator.

Electron-electron interactions effects are ubiquitous in condensed matter systems, especially when charge screening is relatively poor as is the case generally for Dirac-like states in graphene\cite{CastroNeto09,Wehling11} or TIs and TCIs.\cite{Coulomb} In this connection, we discuss possible symmetry breaking orders generated under short-range repulsion $U$ as a first step toward understanding correlation driven effects on Dirac states in the presence of applied strains. We numerically obtain the zero-temperature phase diagram as a function of $U$ and the particle density $n$ and delineate its evolution with increasing strength of strain. Our results not only give insight into the nature of competing orders, but also build the foundation for what may be called ``straintronics'' applications driven by electron correlations.

The paper is organized as follows. In Sec. II, we briefly discuss the effective four-band model for the (001) surface states in the low-energy regime, and its extensions for addressing strain effects in TCIs. Sec. III turns to consider strain effects on properties of the Dirac states, including the interplay between applied strains and perpendicular Zeeman fields. These results then allow us to investigate in Section IV the possible interaction-driven competing orders in Dirac states in the absence as well as presence of strains. Finally, Section V comments on potential applications and implications of our study, and concludes with a summary of our results.

\section{Effective four-band $k\cdot p$ model for the surface states}
We start by reviewing the four-band $k\cdot p$ model for the (001) surface states in a TCI developed previously in Ref.\onlinecite{PhysRevB.87.235317}, and discuss its generalization to account for effects of strains. On the (001) surface, the low-energy surface states can be viewed as two sets of interacting coaxial Dirac cones, originating from the interface between the inverted bands of the TCI and the vacuum: one set is centered at $\bar{X}$ and the other at $\bar{Y}$ in the surface Brillouin zone (SBZ), see Fig.~\ref{fig:lattice} (b). Defining [110], [$\bar{1}$10], and [001] as $x$, $y$, and $z$ directions [see Fig.~\ref{fig:lattice} (a)], respectively, the effective model $H_{\bar{X}}$ around $\bar{X}$ must obey the following three essential symmetries which leave $\bar{X}$ invariant: the mirror reflection about the $xz$-plane, the mirror reflection about the $yz$-plane, and time-reversal symmetry with the corresponding symmetry operations, i.e.
\begin{eqnarray}
M_{xz}H_{\bar{X}}(k_x,k_y)M_{xz}^{-1}=H_{\bar{X}}(k_x,-k_y),\\
M_{yz}H_{\bar{X}}(k_x,k_y)M_{yz}^{-1}=H_{\bar{X}}(-k_x,k_y),\\
TH_{\bar{X}}(k_x,k_y)T^{-1}=H_{\bar{X}}(-k_x,-k_y),
\end{eqnarray}
where $\mathbf{k}=(k_x,k_y)$ is measured from $\bar{X}$.

Informed by the orbital characters of surface states as revealed by first-principles calculations,\cite{PhysRevB.87.235317} one may choose, for instance in the case of SnTe, $\{|p_z,\uparrow;\text{Sn}\rangle$, $
|p_z,\downarrow;\text{Sn}\rangle$, $|p_x,\uparrow;\text{Te}\rangle$, $
|p_x,\downarrow;\text{Te}\rangle\}$ as the basis states of $H_{\bar{X}}$ because one of the interacting Dirac cones is $p_z$-orbital (Sn) rich, while the other is $p_x$-orbital (Te) rich. As a result, the symmetry operations can be represented by the matrices, $M_{xz}=-i\Sigma_{02}$, $M_{yz}=-i\Sigma_{31}$, and $T=-i\Sigma_{02}K$, where $\Sigma_{\alpha\beta}\equiv\sigma_{\alpha}\otimes s_{\beta}$ with the Pauli matrices $\vec{\sigma}$ and $\vec{s}$ acting on orbital and spin spaces, respectively, and $K$ denotes complex conjugation. These considerations restrict $H_{\bar{X}}$ to the form
\begin{eqnarray}
H_{\bar{X}}&=& m\Sigma_{30}+m^\prime\Sigma_{22}+k_x(v_{1x}\Sigma_{02}+v_{2x}
\Sigma_{20}+v_{3x}\Sigma_{32})\nonumber \\
&+& k_y(-v_{1y}\Sigma_{01}+v_{2y}\Sigma_{13}-v_{3y}\Sigma_{31}),\label{H}
\label{eq:surfaceH}
\end{eqnarray}
up to first order in $|\mathbf{k}|$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=8cm]{Fig1_structure.jpg}
\caption{(Color online) (a) Top view of the crystal structure of SnTe. Two types of atoms are drawn in different colors. $x$ ($y$) denotes the direction along [110] ([$\bar{1}$10]). (b) Surface Brillouin zone of SnTe. Without strain, there are four gapless Dirac cones, protected by the mirror symmetry either with respect to $\bar{\Gamma}\bar{X}$ or $\bar{\Gamma}\bar{Y}$.}
\label{fig:lattice}
\end{center}
\end{figure}

For $m^\prime$=0 and $v_{2x}$=$v_{3x}$=$v_{2y}$=$v_{3y}$=0, $H_{\bar{X}}$ becomes block diagonal, reflecting the structure of the two underlying parent Dirac cones with two associated Dirac points, $E_{d1}=-m$ ($p_z$-rich) and $E_{d2}=m$ ($p_x$-rich) at $\bar{X}$ [see Fig.~\ref{fig:b1} (a)]. A non-vanishing value of any of the parameters $m'$, $v_{2x}$, $v_{3x}$, $v_{2y}$ and $v_{3y}$ induces interaction between the two cones, leading to two immediate consequences: (1) The original Dirac points are shifted to $E_{d1(2)}=\mp\sqrt{m^2+m^{\prime 2}}$; and (2) All degeneracies along the intersection of the two cones are lifted except for the two points of a time-reversal symmetric pair on the $\bar{\Gamma}\bar{X}$ line, indicating the emergence of two child Dirac cones protected by the mirror reflection about the $xz$-plane [see Fig.~\ref{fig:lattice} (b)]; hereafter, referred to as the low-energy Dirac points associated with $E_d$=0.

\begin{figure*}[th]
\begin{center}
\includegraphics[width=18cm]{unstrainband.jpg}
\caption{(Color online) (a) Band dispersion of the (001) surface of SnTe along high symmetry lines obtained by setting model parameters in Eq.~(\ref{eq:surfaceH}) as: $a=6.327$ \AA, $m = -0.3$ eV,  $m^\prime= -0.15$ eV, $v_{2x} = v_{2y} = v_{3x} = v_{3y} = 0$ eV\AA, $v_{1x} = -2.3$ eV\AA, and $v_{1y} = -6.5$ eV\AA. $E_d$ and $E_{d1(2)}$ denote energies of the low-energy Dirac points and the two associated parent Dirac points, respectively. $E_{VHS+(-)}$ are energies of the van Hove singularities. (b) and (c) separately show spin textures of valence and conduction bands. Arrows and their lengths indicate directions and magnitudes of the in-plane spin polarization; the color map indicates the distribution of the out-of-the-plane spin polarization. Solid lines mark contours of constant energy.}
\label{fig:b1}
\end{center}
\end{figure*}

Our effective four-band $k\cdot p$ model correctly capture two key features of surface states of TCIs. Firstly, for $E>0$ topology of the constant energy contours changes from two separate Fermi circles at low energies to two concentric ellipses at high energies. This Lifshitz transition\cite{Hsieh12} indicates the presence of van Hove singularities (VHSs) in the underlying electronic spectrum, which are located along the $\bar{X}\bar{M}$ line at
\begin{equation}
E_{VHS+} = |m^{\prime}|,\quad \mathbf{k} = (0,\pm\frac{m}{v_{1y}}),
\end{equation}
if $v_{2x}$, $v_{3x}$, $v_{2y}$, and $v_{3y}$ are neglected. A similar situation arises for $E<0$ with $E_{VHS-}=-E_{VHS+}$. In all cases, we have a logarithmically diverging density of states proportional to $\ln \frac{\Lambda}{|\omega|}$ with $\omega$ and $\Lambda$ representing an energy scale away from the VHS and a cutoff energy for surface states, respectively. Secondly, there is no out-of-the-plane spin component [Fig.~\ref{fig:b1} (b)-(c)] by using symmetry arguments in that the net out-of-the-plane spin polarization, $\langle\Sigma_{03}\rangle=\langle M_{xz}M_{yz}T\Sigma_{03}T^{-1}M_{yz}^{-1}M_{xz}^{-1}\rangle=-\langle \Sigma_{03}\rangle$, yielding  $\langle\Sigma_{03}\rangle =0$.

We turn next to generalize our four-band effective model to include low-energy states throughout the SBZ in the presence of inter-cone interactions. This can be done by noticing that the surface states in the vicinity of $\bar{X}$ and $\bar{Y}$ are related by $C_4$ symmetry, so that the effective model around $\bar{Y}$ can be obtained explicitly by a $C_4$ rotation,
\begin{equation}
H_{\bar{Y}}(k_x,k_y)=\hat{C}_4 H_{\bar{X}}(k_y,-k_x)\hat{C}_4^{-1},
\end{equation}
where $\hat{C}_4=\sigma_0\otimes e^{-i\frac{\pi}{4}s_3}$ and the $p_x$ orbital for Te atoms in the basis states is now replaced by $p_y$ orbital in $H_{\bar{Y}}$. The total surface Hamiltonian then is
\begin{equation}
H_{(001)}=H_{\bar{X}}(k_x,k_y)\oplus H_{\bar{Y}}(k_x,k_y).
\end{equation}
Although $H_{(001)}$ describes an even number of Dirac cones at low energies like the case of a weak TI\cite{Fu07a,Fu07b,LFTCI2,Yan12}, note that the four Dirac points in a TCI do not locate at time-reversal invariant momenta. This key difference can lead to a rather different phase diagram in the presence of the electron-electron interactions\cite{Liu2012906} in a TCI, as discussed in Sec. IV below.

\section{STRAIN EFFECTS}

\subsection{Gap opening and spin textures}
We now discuss effects on the surface states due to strains, which could be either intrinsic or extrinsic. Since the gapless surface states in the Pb$_{1-x}$Sn$_x$(Te,Se) TCIs are mainly protected by mirror symmetries, a perturbation, which breaks one of these symmetries can be expected to at least partly gap the surface spectrum, and modify the spin textures around the associated (massive) Dirac points. Although this is generally true, we will see below that this is not always the case due to other symmetry considerations. For this purpose, we will carry out a systematic analysis of strain induced effects along the lines of Sec. II above.

\begin{figure*}[th]
\begin{center}
\includegraphics[width=18cm]{gap1_v2.jpg}
\caption{(Color online) (a) Band dispersion of (001) surface of SnTe along high symmetry lines with the strain perturbation $\Delta\Sigma_{23}$; parameters are the same as those used in Fig.~\ref{fig:b1} and $\Delta=0.01$ eV. (b) and (c) are spin textures of the conduction band in the presence of the strain and a perpendicular Zeeman field [modeled by adding $H^\prime=h_z\Sigma_{03}$ in Eq.~(\ref{eq:surfaceH}) with $h_z=0.01$ eV], respectively. Arrows and their lengths indicate the direction and the magnitude of in-plane spin polarizations; color maps show the distribution of the out-of-the-plane spin polarization. Solid lines mark constant energy contours.}
\label{fig:f1}
\end{center}
\end{figure*}

A general strain can be described in terms of a symmetric strain tensor $\varepsilon_{ij}=\partial_j u_i$, where $\mathbf{u}$ denotes the displacement field and $i,j=x,y,z$ are chosen to coincide with the principal crystal axes as defined in Sec. II above. Any effective, strain-induced perturbation term corresponding to $\varepsilon_{ij}$ will transform in the same way as the strain tensor under time-reversal and mirror symmetries. By examining the transformation properties of $\Sigma_{\alpha\beta}$ up to first order in $|\mathbf{k}|$, the possible resulting perturbations around $\bar{X}$ are listed in Table \ref{T1}.

\begin{table*}
\caption{Effects of perturbations induced by applied strains. The first row shows the form of each perturbation in terms of a $4\times 4$ matrix in the basis of the effective model around $\bar{X}$. $\Delta$ and $v$ are the coupling strengths. The next two rows indicate if a perturbation breaks certain mirror symmetry. The fourth and fifth rows represent the possible consequences. The last row shows the corresponding strain tensor with which each perturbation can couple.}
\begin{ruledtabular}
\begin{tabular}{ccccccccccccc}
Perturbation &$\Delta\Sigma_{23}$&$vk_x\Sigma_{03}$&$vk_x\Sigma_{11}$&$vk_x\Sigma_{33}$
&$\Delta\Sigma_{10}$&$vk_x\Sigma_{12}$&$\Delta\Sigma_{21}$& $vk_x\Sigma_{01}$&$vk_x\Sigma_{13}$&$vk_x\Sigma_{31}$ \\
Break $xz$-mirror symmetry
&$\checkmark$&$\checkmark$&$\checkmark$&$\checkmark$&$\times$&$\times$
&$\checkmark$&$\checkmark$&$\checkmark$&$\checkmark$ \\
Break $yz$-mirror symmetry
&$\times$&$\times$&$\times$&$\times$&$\checkmark$&$\checkmark$&$\checkmark$
&$\checkmark$&$\checkmark$&$\checkmark$ \\
Open gaps at Dirac points
&$\checkmark$&$\checkmark$&$\checkmark$&$\checkmark$&$\times$&$\times$&$\times$
&$\times$&$\times$&$\times$ \\
Induce out-of-plane spins
&$\checkmark$&$\checkmark$&$\checkmark$&$\checkmark$&$\times$&$\times$&$\times$
&$\times$&$\times$&$\times$ \\
Coupled strain tensor
&$\varepsilon_{yz}$&$\varepsilon_{yz}$&$\varepsilon_{yz}$&$\varepsilon_{yz}$
&$\varepsilon_{xz}$&$\varepsilon_{xz}$&$\varepsilon_{xy}$&$\varepsilon_{xy}$
&$\varepsilon_{xy}$&$\varepsilon_{xy}$
\label{T1}
\end{tabular}
\end{ruledtabular}
\end{table*}

As expected, the strain induced perturbations in the first four columns of Table I can open gaps at both low-energy Dirac points along $\bar{\Gamma}\bar{X}$ due to $M_{xz}$ broken symmetry. More significantly, the resulting spin texture becomes hedgehog-like at low energies, similar to the spin reorientation phenomenon found on the surfaces of manganese-doped Bi$_2$Se$_3$ thin films\cite{Xu12b}.

As a concrete example, let us add the perturbation $\Delta\Sigma_{23}$ in Eq.(\ref{eq:surfaceH}), where $\Delta$ denotes the electron-phonon coupling strength. Clearly, in Fig.~\ref{fig:f1}(a), a gap opens at the Dirac point with its magnitude proportional to $|\Delta|$. The spin-texture of the conduction band of the massive Dirac cone is also depicted in Fig.~\ref{fig:f1}(b). As the absolute value of the eigenenergy $|E|$ becomes smaller ({\it i.e.}, closer to the bottom of the upper cone), the induced out-of-the-plane spin component,
\begin{equation}
\langle\Sigma_{03}\rangle
=\frac{2m'\Delta v_{1x}k_x}{E[E^2-(m^2+m'{}^2+v_{1x}^2k_x^2+v_{1y}^2k_y^2+\Delta^2)]},
\label{eq:sz}
\end{equation}
becomes larger, where we have set $v_{2x}=v_{2y}=v_{3x}=v_{3y}=0$ for simplicity. Several points should be noted here as follows. Firstly, since the strain perturbation does not break TRS, the spin texture for the two massive cones around $\bar{X}$ must form time-reversed partners. This can be inferred from Eq.~(\ref{eq:sz}) by changing the sign of $k_x$. Secondly, the out-of-the-plane spin component of the lower cone is basically opposite to that of the upper cone at the same $|\mathbf{k}|$. Finally, the whole spin texture resulting from the strain differs from the case where the hedgehog-like texture is induced by the perpendicular Zeeman field when both massive Dirac cones around $\bar{X}$ are considered, see Fig.~\ref{fig:f1}(c).

Other strain perturbations in Table \ref{T1}, which break either the $yz$-mirror symmetry or both the mirror symmetries, neither open a gap at the Dirac point nor induce an out-of-the-plane spin component. The robustness of this gapless Dirac point originates from an underlying symmetry: the former type of perturbation is due to the presence of $M_{xz}$ symmetry, while the latter type (shear deformation) is protected by a ``space-time'' symmetry, $C_{2T}=i\Sigma_{31}K$, with a rotation $C_2$ followed by a time-reversal operation. The space-time protection now allows persistence of gapless Dirac points sitting at generic $\mathbf{k}$ points away from the mirror line $\bar{\Gamma}\bar{X}$. This special feature has been noted previously in Refs.\onlinecite{Fang14,Serbyn14}, and it leads to a new type of TCI\cite{Fang15}.

Vanishing out-of-the-plane spin polarization can be proven via symmetry arguments. For instance, consider the surface states around $\bar{X}$ under the strain perturbation $\Delta\Sigma_{10}$. The corresponding Hamiltonian now reads: $\tilde{H}_{\bar{X}} = H_{\bar{X}}+\Delta\Sigma_{10}$. Assuming $|\Psi\rangle$ is an eigenstate of $\tilde{H}_{\bar{X}}$, the $xz$-mirror symmetry guarantees that
\begin{equation}
0=\langle\Psi|[\Sigma_{02},\tilde{H}_{\bar{X}}]|\Psi\rangle=2iv_{1y}k_y
\langle\Psi|\Sigma_{03}|\Psi\rangle.
\end{equation}
Thus, the out-of-the-plane spin component $\langle\Psi|\Sigma_{03}|\Psi\rangle$ vanishes everywhere around $\bar{X}$. Other perturbations can be analyzed in a similar manner.

There other cases deserve some comment. The first two cases involve a uniform expansion and a uniaxial stretch ($C_4$ breaking), which correspond to $\varepsilon_{xx}+\varepsilon_{yy}$ and $\varepsilon_{xx}-\varepsilon_{yy}$, respectively\cite{Zeljkovic15}. Because both these cases respect TRS as well as the mirror symmetries, their net effect is only to renormalize the parameters in the original $H_{\bar{X}}$. Consequently, these perturbations do not open a gap at Dirac points, although positions of the Dirac points could shift in opposite directions along the $\bar{\Gamma}\bar{X}$ line. The third case is an experimentally observed ferroelectric-like distortion\cite{OkadaTCIstrain}, in which two kinds of atoms are displaced along a certain direction in an opposite manner. Denoting the displacement vector $\textbf{d}=(d_x,d_y)$, the non-vanishing component of this distortion, $d_x$ ($d_y$), preserves the mirror symmetry with respect to the principal axis $x$ ($y$), but breaks the rotation and the mirror symmetries along a perpendicular direction. To the zeroth order in $\mathbf{k}$ around $\bar{X}$, the perturbation due to such distortion can be straightforwardly shown to have the form:
\begin{equation}
H_F^\prime = \Delta_{Fx}d_x\Sigma_{10}+\Delta_{Fy}d_y\Sigma_{23},
\end{equation}
where $\Delta_{Fx}$ and $\Delta_{Fy}$ denote the coupling strengths. In fact, the former term is similar to the effect of strain $\varepsilon_{xz}$, preserving the gapless Dirac points, while the latter term is similar to the effect of strain $\varepsilon_{yz}$, resulting in opening a gap instead.

\begin{figure}[tbh]
\begin{center}
\includegraphics[width=7cm]{bandinverse.jpg}
\caption{(Color online) Evolution of band gaps in the Dirac cones under strain perturbation $V_s\Sigma_{23}$ with a fixed Zeeman field, where $h_z=0.1$. The strength of the strain perturbation is given by $V_s=$ (a) 0, (b) 0.03, (c) 0.05, and (d) 0.08. From (b) to (d), one of the Dirac cones around $\bar{Y}$ closes and reopens again due to a band inversion. The model parameters are the same as those used in Fig.~\ref{fig:b1}.}
\label{fig:biv}
\end{center}
\end{figure}

\subsection{Tuning the Chern number via the interplay of strain and a Zeeman fields}
Like a perpendicular Zeeman field, the strain field on the surface of a TCI may not only give rise to spin reorientation, but it may also induce a charge gap in the Dirac cones. However, there are key differences in the effects of strain and Zeeman fields. An applied Zeeman field, which respects $C_2$ rotational symmetry around the $z$-axis, opens gaps with the {\it same} sign in the pair of Dirac cones around $\bar{X}$, while an applied strain field, which breaks $C_2$ rotational symmetry, would induce gaps with {\it opposite} sign. This observation provides the foundation for tenability of the Chern number via the interplay of applied strain and Zeeman fields, suggested in Ref.\onlinecite{Fang14} using a simplified two-band model for each Dirac point.

As a concrete example based on our four-band model, consider an applied Zeeman field along the $z$ direction along with a strain field which breaks $M_{yz}$ but preserves $M_{xz}$ around $\bar{Y}$ on the (001) surface. From symmetry considerations, the strain and Zeeman fields around $\bar{X}$ and $\bar{Y}$ become coupled as
\begin{eqnarray}
H_{\bar{X}} &+& h_z\Sigma_{03}, \nonumber \\
H_{\bar{Y}} &+& h_z\Sigma_{03}+V_s\Sigma_{23},
\end{eqnarray}
where $h_z$ ($V_s$) is the field strength of Zeeman (strain) field. Note that for a given $h_z$, all four low-energy Dirac points open gaps in the absence of an applied strain [see Fig.~\ref{fig:biv}(a)]. When a non-vanishing $V_s\ll \frac{1}{2}h_z$ is introduced, the gaps of the two massive Dirac cones around $\bar{Y}$ evolve in an opposite manner in the sense that one cone increases while the other decreases. This is consistent with a picture in which the strain field induces out-of-the-plane spin components around the two massive Dirac cones which form time-reversed partners with respect to their spin textures [see Fig.~\ref{fig:f1}(b)], and therefore, respond oppositely to the existing $h_z$ field. With increasing strength of $V_s$, the gap in the Dirac cone with a decreasing gap continues to decrease further, becomes gapless at $V_s = \frac{1}{2}h_z$, and reopens inverted again, as shown in Figs.~\ref{fig:biv}(b)-(d). This band inversion indicates that the system undergoes a topological phase transition.

The topological nature of a ground state with broken time-reversal symmetry can be characterized by the Chern number. A non-zero value of the Chern number, $C$, given by
\begin{align}
C=\frac{-1}{\pi }\sum_{n\in\text{occ.}}\sum_{m\neq n}\frac{\text{Im}[\langle u_n|\frac{\partial H}{\partial k_x}|u_m\rangle\langle u_m|\frac{\partial H}{\partial k_y}|u_n\rangle]}{(E_n-E_m)^2},
\end{align}
where the summation over momentum space is implicit, indicates a finite Hall conductance $\sigma=C\frac{e^2}{h}$, which can be obtained by integrating the Berry curvature of the wave functions of the occupied bands over the momentum space\cite{Andrei13}. Fig.~\ref{fig:Chern} presents the Chern number of the system as a function of the Zeeman field strength $h_z$ and the strain strength $V_s$ within the framework of our four-band model. In the presence of a Zeeman field, by varying the strength of strain from small to large values is seen to drive the system from $C=2$ to $C=1$, demonstrating the tunability of the system to the strain field. Note that in our consideration here we have implicitly assumed that the sample is thick enough so that the hybridization between the top and bottom surfaces of the sample can be neglected; with $h_z\neq 0$ the bottom surface will contribute another Chern number $C=2$, which is not shown in Fig.~\ref{fig:Chern} for simplicity.

\begin{figure}[tbh]
\begin{center}
\includegraphics[width=200pt]{Chern_v1.jpg}
\caption{(Color online) Chern number for our four-band model of a TCI as a function of $V_s$ and $h_z$. Model parameters are fixed as: $m$=-2, $m'$=-1, $v_{1x}$=-1, $v_{1y}$=-1.5, with other parameters taken to be zero. The phase boundary (blue line) is given by $h_z=2V_s$.}
\label{fig:Chern}
\end{center}
\end{figure}

\section{EFFECTS OF ELECTRON-ELECTRON INTERACTION}

In discussing electron-electron interactions, we consider short-range repulsive interactions between electrons on the (001) surface of a TCI, and focus on delineating how the phase diagram evolves in the presence of a strain, which breaks one of the mirror symmetries, say, $M_{yz}$, based on the the interacting Hamiltonian:
\begin{equation}
\hat{H}=\sum_{\mathbf{k}}\hat{\psi}^{\dagger}(\mathbf{k})H_{(001)}\hat{\psi}(\mathbf{k})+\hat{H}_{int},
\end{equation}
where
\begin{eqnarray}
\hat{\psi}^{\dagger}(\mathbf{k})&=&[c^{\dagger}_{\bar X,p_z\uparrow}(\mathbf{k}),c^{\dagger}_{\bar X,p_z\downarrow}(\mathbf{k}),
c^{\dagger}_{\bar X,p_x\uparrow}(\mathbf{k}),c^{\dagger}_{\bar X,p_x\downarrow}(\mathbf{k}),\nonumber \\
&& c^{\dagger}_{\bar Y,p_z\uparrow}(\mathbf{k}),c^{\dagger}_{\bar Y,p_z\downarrow}(\mathbf{k}),
c^{\dagger}_{\bar Y,p_y\uparrow}(\mathbf{k}),c^{\dagger}_{\bar Y,p_y\downarrow}(\mathbf{k})],
\label{eq:basis}
\end{eqnarray}
with the subscript $\bar X$ ($\bar Y$) denoting the momentum point involved in the $\mathbf{k}$ expansion. For $H_{(001)}$, we have added the strain perturbation and hence, $H_{(001)}=H_{\bar{X}}\oplus (H_{\bar{Y}}+V_s\Sigma_{23})$. The interaction term, $\hat{H}_{int}=\hat{H}_{U}+\hat{H}_{V}$; $\hat{H}_{U}$ and $\hat{H}_{V}$ denote repulsive contact interactions between like and unlike orbitals, respectively, with
\begin{eqnarray}
\hat{H}_{U} &=& U\int d^2 r\;\sum_{\eta=p_z,p_x,p_y}\hat{n}_{\eta\uparrow}({\mathbf{r}})
\hat{n}_{\eta\downarrow}({\mathbf{r}}), \nonumber \\
\hat{H}_{V} &=& \frac{V}{2}\int d^2 r\;\sum_{\eta\neq\eta^\prime}\hat{n}_{\eta}({\mathbf{r}})
\hat{n}_{\eta^\prime}({\mathbf{r}}),
\label{Hint}
\end{eqnarray}
where $\hat{n}_{\eta}(\mathbf{r})=\sum_{s}\hat{n}_{\eta s}(\mathbf{r})=\sum_{\mathbf{Q},s}\tilde{c}^\dagger_{\mathbf{Q},\eta s}\tilde{c}_{\mathbf{Q},\eta s}(\mathbf{r})$ ($\mathbf{Q}=\bar{X},\bar{Y}$, $\eta=p_x,p_y,p_z$) [$\tilde{c}_{\mathbf{Q},\eta s}(\mathbf{r})$ are defined in Eq.~(\ref{eq:FTc}) below].

Since we are treating only the low-energy, long-wavelength physics, it is reasonable to consider field operators $c_{\mathbf{Q},\eta\sigma}(\mathbf{r})$, which vary slowly on the scale of the lattice constant. This can be done by Fourier transforming the operators:
\begin{equation}
\tilde{c}_{\mathbf{Q},\eta s}(\mathbf{r}) = \sum_{\mathbf{K}}e^{i\mathbf{K}\cdot\mathbf{r}}
c_{\mathbf{Q},\eta s}(\mathbf{K}) = e^{i\mathbf{Q}\cdot\mathbf{r}}
c_{\mathbf{Q},\eta s}(\mathbf{r}),
\label{eq:FTc}
\end{equation}
where $\mathbf{K}=\mathbf{Q}+\mathbf{k}$, expanding from the origin (0,0). In terms of these field and density operators, and the identities given in the Appendix, the full Hamiltonian $\hat{H}$ can be rewritten as
\begin{equation}
\hat{H}=\int d^2 r\;\hat{\Psi}^{\dag}H_{(001)}\hat{\Psi}+\hat{H}_{int},
\end{equation}
where $\hat{\Psi}^{\dagger}(\mathbf{r})$ has the same form as Eq.~(\ref{eq:basis}) except that the $\mathbf{k}$-dependence is now replaced by $\mathbf{r}$.

In order to determine the ground state of the system as a function of the interaction strength and fermion density, we employ the self-consistent mean-field (MF) approach, which is expected to be reliable as long as the interaction strength is much smaller than the bulk band gap. Taking $U=V$ as a representative case, we decouple the $\hat{H}_{int}$ into bilinear fermion terms. After some straightforward but tedious algebra, we obtain 52 order parameters as well as 12 renormalized band parameters as detailed in the Appendix, which can be generally expressed as
\begin{equation}
O_{\gamma\alpha\beta}=\langle\hat{\Psi}^{\dagger}\Sigma_{\gamma\alpha\beta}
\hat{\Psi}\rangle
\label{eq:order}
\end{equation}
where a summation over momentum and occupied states is implicit, and $\Sigma_{\gamma\alpha\beta}=\tau_{\gamma}\otimes\Sigma_{\alpha\beta}$, with the Pauli matrices, $\tau_{\gamma}$, acting on $\bar{X}\bar{Y}$ pseudo-spin (valley) space. Note that, here we consider only the symmetry breaking orders in the particle-hole channel: superconductivity in the particle-particle channel due to weak onsite repulsion can only be achieved via beyond MF treatment, usually with exponentially small transition temperature\cite{Kohn65}.

In this connection, we define the ``particle density'' at a given chemical potential $\mu$ as the density deviation from the total particle density where $\mu=E_d=0$, namely,
\begin{equation}
n=\int_{k<k_{\Lambda}}^{E(\mathbf{k})=\mu} \frac{d^2k}{(2\pi)^2}-\int_{k<k_{\Lambda}}^{E_0(\mathbf{k})=0} \frac{d^2k}{(2\pi)^2},
\label{eq:density}
\end{equation}
where $E_0(\mathbf{k})$ ($E(\mathbf{k})$) is the eigenenergy of the non-interacting (MF-decoupled interacting) system; $k_{\Lambda}$ is a cutoff momentum, which is chosen such that our results are insensitive to its magnitude. Note that by any induced valley ``polarization'', $O_{300}$, we mean
\begin{equation}
O_{300}=\langle\hat{\Psi}^{\dagger}\Sigma_{300}
\hat{\Psi}\rangle
-\int^{E_0(\mathbf{k})=0}_{k<k_{\Lambda}}\frac{d^2k}{(2\pi)^2}\langle\Sigma_{300}\rangle.
\label{eq:orbital}
\end{equation}
That is, the ``valley polarization'' due to external strain in the absence of interactions is subtracted in defining $O_{300}$.

\begin{figure}
\begin{center}
\includegraphics[width=200pt]{phasediagram_1.jpg}
\caption{(Color online) $U$ versus $n$ phase diagram without strain. $n_{\text{VHS}}$ denotes the particle density at $E_{VHS+}$. Solid and dashed lines mark phase boundaries of continuous and abrupt phase transitions, respectively.}
\label{fig:PD1}
\end{center}
\end{figure}

For a given $U$ and $n$, the ground state can then be determined by minimizing the MF free energy with respect to various order parameters using the form detailed in the Appendix. We numerically solve the resulting set of coupled equations self-consistently to obtain the zero-temperature $U$-$n$ phase diagrams both with and without the applied strain. Note that, in order to make the underlying physics more transparent, we use the following model parameters in computations: $m=-2$, $m'=-1$, $v_{1x}=-1$, $v_{1y}=-1.5$, and a cutoff energy, $\Lambda\equiv E_{0}(k_{\Lambda})=2\sqrt{m^2+m'{}^2}$, which determines $k_{\Lambda}$. However, we expect our results to be generic and relevant more generally, detailed effects of material-specific parameters notwithstanding.

\subsection{System without strain}
We begin by considering a system without strain, {\it i.e.}, $V_s=0$. In the order parameter space we explored in this case, the leading orders (aside from the metal phase) are basically the ferromagnetic state (FM), $O_{003}$, and the $p_z$-orbital rich $(\pi,\pi)$ charge density wave ($p_z$-CDW):
\begin{equation}
\frac{O_{100}+O_{130}}{2}=\Big\langle\sum_{\mathbf{r}}(-1)^{x+y}[\hat n_{p_z,\uparrow}(\mathbf{r})+\hat n_{p_z,\downarrow}(\mathbf{r})]\Big\rangle,
\end{equation}
where the length scale is in units of the lattice constant.  The zero-temperature $U$-$n$ phase diagram is shown in Fig.~\ref{fig:PD1}. The $p_z$-CDW phase corresponds to a non-vanishing $\frac{O_{100}+O_{130}}{2}$ (the leading piece), and the FM phase represents a non-vanishing value of $O_{003}$ (the leading piece). Finally, in the metal phase, all symmetries are preserved with vanishing values of all 52 order parameters.

At $n=0$, the ``$p_z$-CDW'' phase appears when the interaction strength $U$ is larger than the critical value $U_c\approx 2.7$. This phase involves a finite value of $\frac{O_{100}+O_{130}}{2}$; it is associated with broken translation symmetry, while the TRS and both the mirror symmetries remain intact [see schematic Fig.~\ref{fig:realspace}(b)], so that each Dirac cone remains gapless. The existence of a critical value of $U$ can be understood by noticing that in the non-interacting limit, there are only four Dirac points with zero density of states at the Fermi level. The necessity of a threshold value of $U$ has also been predicted theoretically in 2D systems with linear energy dispersion, such as graphene and the surfaces of 3D weak TIs\cite{Raghu08,Liu2012906}.

Upon electron doping, the system assumes either the FM phase [see Fig.~\ref{fig:realspace}(a)] at relatively high particle densities, or the $p_z$-CDW phase [see Fig.~\ref{fig:realspace}(b)] at low particle densities. The FM phase is associated with broken TRS as well as both the mirror symmetries simultaneously, resulting in gaps to open up at all the Dirac points, although the corresponding ground state is still metallic with finite $\mu$. Note that the critical $U$, above which the system enters into the gapless $p_z$-CDW phase, is less than the critical value for the FM phase, indicating a favorable free energy gain compared to the competing FM phase at low doping.

\begin{figure}[th]
\begin{center}
\includegraphics[width=8cm]{schematic_phase.jpg}
\caption{(Color online) Schematic plots illustrating (a) the ferromagnetic state (FM), and (b) the $p_z$-orbital rich, $(\pi,\pi)$ charge density wave state ($p_z$-CDW). Different colors label the two distinct atoms involved; arrows depict the net out-of-the-plane spin polarization.}
\label{fig:realspace}
\end{center}
\end{figure}

We further note several points in connection with the phase diagram of Fig.~\ref{fig:PD1} as follows: (1) The spin density wave (SDW) phase, as anticipated usually in Dirac systems like graphene, does not occur here at $n=0$ for large $U$. This could be attributed to the non-trivial orbital textures and strong spin-orbit coupling in our case; (2) We find that the dynamically generated strain-like order, which could gap the spectrum, is relatively disfavored as the system tries to preserve the combined $C_{2T}$ symmetry, although such a phase becomes favored in the presence of external strain; (3) If we suppress inter-valley scattering, our preliminary results suggest that the CDW phase is also suppressed\cite{long-range}; and finally, (4) As we pointed out already in Sec. II above, the density of states at an energy scale $\omega$ near a van Hove singularity diverges like $\ln |\omega|$, indicating propensity for the occurrence of more symmetry breaking orders (including those in the particle-particle channel). The present mean-field treatment would then be inappropriate, even in the weak-coupling limit.

\subsection{System under moderate strain}
Here we consider an applied strain with a moderate coupling strength, $V_s=|m'|$ [of order $\mathcal{O}(E_{VHS+}$)], which breaks the $yz$-mirror symmetry. The resulting zero-temperature phase diagram is shown in Fig.~\ref{fig:PD2}. In contrast to the case without strain, a coexisting phase lying between the $p_z$-CDW and FM phases is now seen to emerge. Suppression of the $p_z$-CDW phase with strain can be anticipated on physical grounds. The reason is that any $(\pi,\pi)$ order requires good nesting between the valleys, but the strain weakens the nesting tendency by opening gaps at Dirac points along the $yz$-mirror line. Notably, in the strained system, spin polarization in the FM phase is mainly contributed by the $p_x$ orbital due to the broken mirror symmetry, which lifts the $p_y$ orbital to higher energies.

\begin{figure}
\begin{center}
\includegraphics[width=200pt]{phasediagram_2.jpg}
\caption{(Color online) $U$ versus $n$ phase diagram in the presence of a moderate strain. Lines mark boundaries between various continuous phases.}
\label{fig:PD2}
\end{center}
\end{figure}

\subsection{System with large strain}
When the coupling strength of the applied strain $V_s$ becomes the largest energy scale in the system, all other degrees of freedom can be ignored at the $\bar Y$ valley as states around $\bar Y$ are all gapped out, leading to the phase diagram of Fig.~\ref{fig:PD3}. All $(\pi,\pi)$ orders are seen to disappear, and only the FM phase survives for $U>U_c(n)$ at a given density $n$. The presence of the concave boundary region separating the metallic and FM phases in Fig. 8 can be understood from energetics: At any given density in this region, the chemical potential lies within the gap induced by the FM order at $E_{d2}$, yielding a gain in free energy and thus lowers $U_c$. As the particle density increases, the chemical potential is eventually unable to locate within the gap and this argument fails. Interestingly, at $n=0$, although the ground state is insulating with a small gap due to TRS breaking, the system does not exhibit quantum anomalous Hall effect. This should be contrasted sharply with the case of a Zeeman field applied perpendicular to the (001) surface of a TCI in the non-interacting limit, where such a broken TRS leads to non-vanishing Chern number as shown in Sec. IIIB above.

\begin{figure}
\begin{center}
\includegraphics[width=200pt]{phasediagram_3.jpg}
\caption{(Color online) $U$ versus $n$ phase diagram in the case of large strain. Solid line marks the boundary between the two continuous phases.}
\label{fig:PD3}
\end{center}
\end{figure}

\section{Discussion and conclusion}

Some of the strains discussed in this study could be realized experimentally through chemical doping. Since PbSe and SnSe assume cubic and orthorhombic structures\cite{structure}, respectively, strains or lattice distortions could likely be manipulated by adjusting chemical compositions. In fact, a ferroelectric-like lattice distortion has been observed in the TCIs via STM/STS spectroscopy\cite{OkadaTCIstrain}. Mechanical tuning could provide another approach. For example, a piezoelectric layer could be grown on the surface of a TCI to allow manipulation of strains via electric field as proposed in Ref.\onlinecite{Fang14}.

As to practical applications, we have shown clearly that the Chern number of the surface states can be tuned through the interplay of applied strain and Zeeman fields. This provides a potential new handle for controlling the topological conducting edge channels for low power-consuming, next generation electronic technologies. Furthermore, our analysis shows that in the presence of electron-electron interactions, ferromagnetism could be turned on/off via strain fields, which, offering a novel pathway toward spintronics applications.

In summary, we have systematically investigated the effects of strains, which break either one or both mirror symmetries without breaking time-reversal symmetry, on the (001) surface of the TCIs. Under mirror-symmetry breaking strains, the low energy Dirac cones not only become massive, but also develop hedgehog-like spin textures. We show that the Chern number of the surface states can be tuned via applied strain fields in the presence of a Zeeman field. Finally, we delineate the competing orders that can result from varying strengths of short-range (repulsive) electron-electron interactions at the mean-field level, and show that the resulting correlated phases are amenable to control through strain fields. The tunability of the interacting as well as the non-interacting electronic structures and topological states of the TCIs via strain fields and their interplay with Zeeman fields as revealed by our study suggests new pathways for developing spintronics and other applications platforms based on the TCIs.

\begin{acknowledgments}
We thank Chen Fang and Hong Yao for useful discussions. W.F.T. and C.Y.H. acknowledge the support from MOST in Taiwan under Grant No.103-2112-M-110-008-MY3. The work at Northeastern University was supported by the US Department of Energy (DOE), Office of Science, Basic Energy Sciences grant number DE-FG02-07ER46352, and benefited from Northeastern University's Advanced Scientific Computation Center (ASCC) and the NERSC supercomputing center through DOE grant number DE-AC02-05CH11231. H.L. acknowledges the Singapore National Research Foundation for the support under NRF Award No. NRF-NRFF2013-03.
\end{acknowledgments}