\section{}
\section{Introduction}
In the past two decades, circuit quantum electrodynamics (circuit QED) has become a field of growing interest for quantum information processing and also to realize new regimes in quantum optics \cite{Roy2017, Gu2017,GonzCira17_2,PaulGonz18,QuijJoha18,Wallraff2004,Blais2004,Wendin2017,Gu2017,Kockum2019}. The restriction to one dimensional (1D) waveguides  in circuit and waveguide QED enhances directionality and reduces losses and therefore has a great advantage over higher dimensional systems to reach strong- and ultrastrong-coupling regimes \cite{ZhenBara10,Tufarelli2013,Tufarelli2014,SancMart14,SancZuec16,SancGarc17,DevoScho07,BambOgaw14,USC_Anton,ZuecGarc19,MlynWall14}.
A typical circuit QED set-up consists of a superconducting qubit coupled to a 1D transmission line (TL) \cite{Wallraff2004,Blais2004,Wendin2017,Gu2017,Pozar}. Superconducting qubits are artificial atoms built with a non-linear Josephson Junction (JJ), that creates an anharmonic energy spectrum \cite{Kockum2019}. There are different kinds of superconducting qubits like flux qubits, phase qubits, and charge qubits \cite{Wendin2017,Clarke2008}.
A 1D transmission line can be modelled by coupled LC oscillators, each having a characteristic impedance of $Z_0 = \sqrt{L_0/C_0} \approx 50\,\si{\Omega}$, smaller than the quantum resistance $R_Q = \hbar/(2e)^2 \approx 1.0 \, \si{k \Omega}$. But recent studies show that it is possible to reach impedances comparable to the quantum resistance or higher \cite{Pechenezhskiy2020,Niepce2019,Masluk2012,Weissl2015,Samkharadze2016,Kuzmin2019,PuertasMartinez2019}. This can be realized by building circuits made of arrays of JJs \cite{Bell2012,Winkel2020,Krupko2018,Masluk2012,Weissl2015,Kuzmin2019,PuertasMartinez2019} or high-kinetic-inductances materials, called superinductors \cite{Peltonen2018,Grunhaupt2018,Samkharadze2016,Niepce2019,Niepce2020,Pechenezhskiy2019}. High impedance JJ arrays and superinductors are for example used in the Fluxonium qubit \cite{Manucharyan2009,Gruenhaupt2019,Nguyen2019,Koch2009}, which has reduced charge noise sensitivity and can have relaxation times up to milliseconds \cite{Pop2014,Nguyen2019}.
This also has an advantage for metrology, since the charge noise insensitivity makes it possible to measure the current very accurately \cite{Guichard2010}.
Furthermore, high-impedance resonators make it possible for light-matter interaction to reach strong coupling regimes, due to strong coupling to vacuum fluctuations \cite{Stockklauser2017}.\\
In this article, we investigate the spontaneous emission of a transmon \cite{Koch2007} capacitively coupled to a 1D TL that is shorted at one end. This system is known as an "atom in front of a mirror" \cite{Hoi2015,Wen2018,Wen2019,Peng2016,FornDiaz2017}. Instead of using a Markovian master equation approach, we are taking the photon traveling time fully into account, making the dynamics non-trivial
\cite{Dorner2002,Guo2017,Tufarelli2014,Bradford2013,Grimsmo2015,Pichler2016,Pichler2017,Guimond_2017,SchnBusc16,GonzCira17,Wiegand2020}. Furthermore, we explore the above-mentioned regime of a TL impedance exceeding the quantum resistance $Z_0 \gg R_Q$. We find that the system behaves qualitatively different, compared to the well studied low impedance regime $Z_0 \ll R_Q$. The atom reflects strongly at all frequencies, except its transition frequency.  Together with the mirror, it thus forms a cavity and when the transition frequency is close to a cavity mode, we find a vacuum Rabi-splitting, resulting in Rabi oscillations in the spontaneous emission. In this regime, all dynamical timescales are independent of the coupling capacitance and instead depend on the intrinsic transmon capacitance and the TL impedance.
Another cavity-free system that shows Rabi splitting is an artificial atom coupled surface acoustic waves \cite{Ask2019}.
\section{Circuit-QED model}
\label{circuitQED}

\begin{figure}[t]
	\centering
	\begin{minipage}{1\linewidth}
		\begin{overpic}[width=1\linewidth]{TL2.pdf}
			\put(51,20){ $C_c$}
			\put(46,17){ $\phi_J$}
			\put(56,27){ $\phi_0$}
			\put(45,10){ $E_J$}
			\put(69,10){ $C_J$}
			\put(30,23){ $Z_0$}
			\put(78,23){ $Z_0$}
			\put(-3,23){ $a)$}
			\put(-3,-6){ $b)$}

		\end{overpic}
	\end{minipage}
	\begin{minipage}{1\linewidth}
		\vspace{0.2cm}
		\begin{overpic}[width=1\linewidth]{Sketch2OneColour.pdf}
			\put(75,21){$\VRout$}
			\put(30,21){$\VLin$}
			\put(75,5){$\VRin$}
			\put(30,5){$\VLout$}
			\put(75,0){$L$}
		\end{overpic}
	\end{minipage}

	\caption{a) The circuit model of a transmon coupled through the coupling capacitance $C_c$ to a semi-infinite 1D TL with impedance $Z_0$.  The  Josephson energy, flux, and capacitance of the transmon are denoted by $E_J$, $\phi_J$, and $C_J$. The flux on the coupling capacitance $C_c$ is denoted by $\phi_0$ with the corresponding voltage $V_0 = \dot \phi_0$.\\
		b) A sketch of the system depicting an atom in front of a mirror coupled to incoming/outgoing microwave fields to the left/right characterised by their respective voltages $V_{L/R}^{\text{in/out}}$ at the transmon. The mirror couples the fields to the right $\VRin(t)=-\VRout(t-T)$, introducing the time of propagation to the mirror and back $T=2L/v$.
	}
	\label{fig:System}
\end{figure}
Our system consists of a transmon qubit capacitively coupled to a semi-infinite 1D TL at a distance $L$ from its grounded end (see Fig.~\ref{fig:System}). The transmon qubit consists of a superconducting anharmonic $LC$-oscillator, where the inductive ($L$) element is formed by a Josephson junction (JJ) with characteristic energy $E_J$ in parallel with a capacitor ($C$) with capacitance $C_J$. The sinusoidal current-phase relation of the JJ makes the energy spectrum of the transmon qubit anharmonic, allowing for excitation with a single microwave photon using standard harmonic microwave sources. The transmon is capacitively coupled to a microwave TL, characterised by its inductance per unit length $L_0$ and capacitance per unit length $C_0$. The velocity and impedance of the electromagnetic field inside the TL is given by $v= 1/\sqrt{L_0 C_0}$ and $Z_0=\sqrt{L_0/C_0}$, respectively.
Using the standard circuit quantization procedure \cite{Devoret1995,Wiegand2020}, we can derive
the Heisenberg equations of motion for the charge $p_0(t)$ on the coupling capacitor $C_c$, the charge $p_J(t)$ on $C_J$ and its conjugate flux $\phi_J (t)$, giving the phase difference $2e \phi_J (t)/\hbar$ over the JJ. Denoting the operators for the voltages of the incoming/outgoing microwave fields to the left/right of the transmon as $V_{L/R}^{\text{in/out}}(t)$ (see Fig.~\ref{fig:System}b) these equations are

\begin{align}
\partial_t \phi_J (t) &= \frac{1}{C_J} (p_J(t)  + p_0(t) ), \label{eq:EoMphiJa}\\
	\partial_t p_J (t) &= - E_J \frac{2 e}{\hbar} \sin ( \frac{2 e}{\hbar}\phi_J(t) ) \approx -\frac{\phi_J(t)}{L_J}, \label{eq:EoMpJa} \\
	\partial_t p_0(t)  &= \frac{2 p_0(t)}{Z_0 C_\Sigma}  + \frac{2 p_J(t)}{Z_0 C_J}  - \frac{2}{Z_0} \left[ V_L^{\text{in}}(t)+V_R^{\text{in}}(t)\right]   , \label{eq:EoMp0OpenTL}\\
	V_{L/R}^{\text{out}}(t) &= V_{R/L}^{\text{in}}(t) - \frac{Z_0}{2} \partial_t p_0(t), \label{eq:V}
\end{align}
where we denoted the capacitance to ground seen by the JJ as $C_\Sigma = C_c C_J / (C_c + C_J)$ and in the second equation, we introduced the Josephson inductance $L_J =\hbar^2/4e^2E_J$, which describes the linearized dynamics of the Josephson junction. This approximation is obviously good in the weak excitation regime $|\phi_J(t)| < \hbar/2e$ and will also be sufficient to describe the spontaneous emission, where the transmon is initially excited by a single microwave photon

\section{Reflection}

\subsection{Open TL}
To characterise how the behaviour of the system changes when we increase the TL impedance $Z_0$, we first investigate the reflection of microwaves from the transmon coupled to an open TL, i.e. without a mirror.
Since the equations of motion are linear, we can express the reflected field operator $\VLout$ in terms of the incoming probe field operator $\VLin$ by Fourier transforming the equations of motion \eqref{eq:EoMphiJa}-\eqref{eq:V}, assuming no incoming field from the right ($\VRin=0$). The expression for the frequency dependent reflection coefficient is given by

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{Reflection.pdf}
	\caption{Reflection of a transmon in an open TL for different ratios of the TL and qubit impedance $Z_0/Z_J$. The curves show the reflection for $\frac{C_c}{C_c+C_J} = 0.1$ and  $Z_0/Z_J = 0.1$ (purple), $Z_0/Z_J = 1$ (blue), $Z_0/Z_J = 10$ (green), $Z_0/Z_J = 100$ (yellow), and $Z_0/Z_J = 1000$ (red) from low to high  TL impedance. For low impedance, the qubit reflects only at the coupled qubit frequency $\omega_0$, but for high impedance, it reflects everywhere but the uncoupled qubit frequency $\omega_J$.
	}
	\label{fig:Reflection}
\end{figure}
\begin{align}
r(\omega) \equiv \frac{\VLout(\omega)}{\VLin(\omega)} = \frac{C_c Z_0 \omega \left( \frac{\omega^2}{\omega_J^2} - 1 \right) }{  2 i \left( 1 - \frac{\omega^2}{\omega_0^2}\right) + C_c Z_0 \omega \left( \frac{\omega^2}{\omega_J^2} - 1 \right)}
\label{eq:r_open_line}
\end{align}
where $\omega_0 = 1/\sqrt{L_J (C_c + C_J)}$ is the resonance frequency of the coupled transmon and $\omega_J = 1/\sqrt{L_J C_J}$ is the resonance frequency of the bare (uncoupled) transmon. In Fig.~\ref{fig:Reflection} the reflection around the transmon resonance frequencies is shown for different values of $Z_0$. We see that for low impedance $Z_0 C_c \omega < 1$, the reflection is weak except at $\omega_0$ where it is unity, due to resonant reflection from the transmon \cite{Lindkvist2013}. For high impedance $Z_0 C_c \omega > 1$ we instead see strong reflection at all frequencies, except around the "new" resonance frequency $\omega_J$,  where we find zero reflection independent of $Z_0$. The crossover occurs at $Z_0 \sim Z_J C_J/C_c$, introducing the transmon impedance
\begin{align}
Z_J = \sqrt{\frac{L_J}{C_J}} = R_Q \sqrt{\frac{2 E_C}{E_J}},
\end{align}
where $E_C = e^2/(2C_J)$ is the charging energy of the transmon.

In the high impedance regime, the strong scattering away from $\omega_J$ occurs due to the comparably strong capacitive coupling to ground at the transmon {\em without exciting the transmon}. Close to $\omega_J$, the resonantly excited transmon counteracts this capacitive coupling and effectively acts like an open circuit. This is the opposite behaviour compared to the low impedance regime, where the transmon is effectively an open circuit at all frequencies, except at its resonance frequency $\omega_0$, where it acts like a shorted circuit, giving full reflection. By fitting Lorentzians to Eq.~\eqref{eq:r_open_line}, we can extract the high impedance coupling strength $\gamma_J=2/Z_0 C_J$, which in contrast to the low-impedance expression $\gamma_0 = Z_0 C_c^2 /2 L_J (C_J+C_c)^2$ does not depend on $C_c$ and decreases with increasing $Z_0$. In this regime, the voltage $V_0=\VLin+\VLout=\VRin+\VRout$ at the node coupling to the TL oscillates with the full voltage across the JJ, $V_J=\dot{\phi_J}$. Due to the large $Z_0$, the amount of current passing through the TL is small ($\propto V_J/Z_0$), thus not being able to change the voltage $V_0-V_J$ across $C_c$ significantly, i.e. $|V_0-V_J| \ll V_J$. In the low impedance regime, we instead have large currents flowing through the TL, keeping the voltage at the coupling node close to zero, i.e. $|V_0| \ll |V_J|$. These currents obviously scale with $C_c$ and the energy dissipation scales with $Z_0$.
In the following, we investigate how the mirror affects the scattering.
\subsection{Mirror}
\begin{figure}[t]
	\centering

		\begin{overpic}[width=1\linewidth]{FieldPositionLabel3.pdf}
			\put(37,63.7){$\frac{\omega_c}{\omega_J} = 0.98$}
			\put(37,57){$\frac{\omega_c}{\omega_J} = 0.99$}
			\put(37,50.3){$\frac{\omega_c}{\omega_J} = 1$}
			\put(37,43.7){$\frac{\omega_c}{\omega_J} = 1.01$}
			\put(37,36.9){$\frac{\omega_c}{\omega_J} = 1.02$}
			\put(36,32){$\xleftarrow{\hspace*{1.1cm}}$}
	        \put(56.5,32){$\xrightarrow{\hspace*{1.1cm}}$}
	        \put(52.5,32){$\Omega$}
		\end{overpic}
	\caption{  Amplitude of the electromagnetic field between the qubit and the mirror $|f(\omega)| =  \left| \VRout(\omega)/\VLin(\omega) \right|$ for $\frac{C_c}{C_J} = 0.1$ and  $Z_0/Z_J = 1000$. The delay time is given by $T = 2 \pi n /\omega_c, n=1$ and $\omega_c$ varies with respect to $\omega_J$. The different colors show $\omega_c/\omega_J = 0.98$ (red), $\omega_c/\omega_J = 0.99$ (yellow), $\omega_c/\omega_J = 1$ (green), $\omega_c/\omega_J = 1.01$ (blue), $\omega_c/\omega_J = 1.02$ (purple). All curves, show a splitting which can be understood as a Rabi splitting. Without detuning between $\omega_J$ and $\omega_c$, the splitting is centered exactly around $\omega/\omega_J = 1$. The inset shows the appearance of higher cavity modes appearing at $\omega \approx n \omega_c$, where  $\omega_c T= 2 \pi$ for $\omega_c/\omega_J = 1$ corresponding to the green curve in the main figure.
	}
	\label{fig:ReflectionMirror}
\end{figure}
The mirror couples the fields to the right of the transmon $\VRin(t)=-\VRout(t-T)$, introducing the time of propagation to the mirror and back $T=2L/v$, see Fig~\ref{fig:System} b). Similarly as before, we can find the response to a harmonic field incoming from the left by Fourier transformation of the equations of motion.
Since the absolute value of the reflection for the transmon in front of a mirror is always unity, we are now interested in the frequency dependence of the ratio between the trapped field (between the qubit and the mirror) and the incoming field, which is given by $f(\omega) \equiv \VRout(\omega)/\VLin(\omega)$,
\begin{align}
\label{f_omega}
	f(\omega)= \frac{ \left( \frac{\omega^2}{\omega_0^2} - 1 \right) }{    \left( 1 - \frac{\omega^2}{\omega_0^2}\right) - i \frac{C_c Z_0 \omega }{2} \left( \frac{\omega^2}{\omega_J^2} - 1 \right) \left( e^{i \omega T} - 1  \right) },
\end{align}
which is shown in Fig.~\ref{fig:ReflectionMirror}.
In the high impedance regime, we now find cavity resonances between the highly reflective atom and the mirror when the frequency is close to $n\ \omega_c$ for $n=1,2,\dots$ and $\omega_c = 2 \pi/T$, as shown by the peaks in the inset of Fig.~\ref{fig:ReflectionMirror}. These are broadened by the coupling to the TL by $\gamma_c^n=|t(n\ \omega_c)|^2/T$, where $t(\omega)$ is the transmission across the transmon and $|t(\omega)|^2=1-|r(\omega)|^2$.
We find that the effect on the transmon resonance close to $\omega_J$ is simply to reduce its broadening with a factor of two to $\gamma_J^m = 1/Z_0 C_J$, away from any qubit-cavity resonance $\omega_J \approx n\ \omega_c$.
As shown in the main panel of Fig.~\ref{fig:ReflectionMirror}, on resonance $\omega_J \approx n\ \omega_c$ we find an avoided crossing with the coupling strength
\begin{align}
    \Omega &= \frac{2}{ \sqrt{T C_J Z_0 } }
 =\frac{2\omega_J}{ \sqrt{2 \pi n C_J Z_0 \omega_J} }
    = \frac{ 2 \omega_J}{ \sqrt{2 \pi n  \frac{Z_0}{Z_J} } }.
\label{eq:omegaB}
\end{align}
As we will see in the next section, where we investigate the spontaneous emission dynamics of the transmon, we find that this coupling indeed gives rise to vacuum Rabi-oscillations between the transmon and the cavity mode.

In the low impedance regime $Z_0 C_c \omega < 1$, $f(\omega)$ is instead close to unity, indicating only little scattering from the transmon for all frequencies far from the transmon resonance $\omega = \omega_0$. Here, $f(\omega_0)=0$, since the field is reflected by the transmon and does not reach the mirror. When the transmon is located at a distance corresponding to a node of the electromagnetic field at its resonance frequency, i.e. $\omega_0 T = 2 n \pi$, it is in a dark state and is thus completely invisible to the incoming field at frequency $\omega_0$, giving instead $f(\omega_0)=1$. In the dark state, both the transmon and the field between the transmon and the mirror are excited. Thus, if the distance is slightly longer/shorter than the node, the state is no longer completely dark, and we instead get a pronounced scattering resonance ($|f(\omega)|\gg 1$) at frequencies slightly lower/higher than $\omega_0$, see e.g. the purple line in Fig.~\ref{fig:ReflectionMirror}. For higher $Z_0$, we see that this dark state resonance moves in frequency towards the cavity frequency $\omega_c = 2 \pi/T$. As shown in the supplemental material, for small $C_c/C_J \ll 1$, and $\omega_0=n\ \omega_c$ we can find vacuum Rabi-oscillations damped towards a finite dark state population.

 In the following, we investigate how the high-impedance TL influences spontaneous emission of the transmon.
\section{Spontaneous emission and Rabi oscillations}
\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{EnergiesInsetScaled.pdf}
	\caption{Energy of the transmon (blue), Field between the mirror and transmon (pink), outgoing field to the left (orange), sum of all the energies (red), approximated decay $\gamma_J^m$ (cyan, dashed) and approximated Rabi frequency $\Omega$ and decay $e^{-\gamma^m_J t/2} \cos^2 (\Omega t/2)$ (yellow, dashed) for $C_c/C_J = 0.1, Z_0/Z_J = 1000, T = 2 \pi/\omega_C, \omega_C = \omega_J$ as a function of periods of the Rabi oscillations. The Inset shows the energy of the transmon (blue) and the flux $\phi_J^2$ (green) for the first oscillation of the frequency $\Omega$. This clearly shows the different timescales of the energy compared to the variables of the transmon.
	}
\label{fig:Energies}
\end{figure}
We consider the case of a transmon initially excited at time $t=0$ with a finite flux $\phi_J(0)>0$, while the other qubit variables are zero $p_J(0)=p_0(0)=0$ and the TL is in the vacuum state. The qubit energy
\begin{equation}
E_q(t)=\frac{\left(p_J(t)+p_0(t)\right)^2}{2C_J}+\frac{p_0(t)^2}{2C_c} + \frac{\phi_J(t)^2}{2L_J},
\end{equation}
is the sum of the capacitive energy on the two capacitances and the inductive energy in the JJ. The current amplitude emitted from the transmon into the TL is $\partial_t p_0(t)$ and from this we can write the change of the energy $E_R(t)$ of the field between the transmon and the mirror as
\begin{equation}
\partial_t E_R(t) = \frac{Z_0}{4} \left[ (\partial_t p_0(t))^2 - (\partial_t p_0(t-T))^2 \right],
\end{equation}
where the first term corresponds to the instantaneous power emitted into the TL and the second term is the instantaneous power coming back from the mirror. The change of the energy of the field to the left of the transmon $E_L(t)$ is given by the instantaneous left-moving power leaving the system,
\begin{equation}
\partial_t E_L(t) = \frac{Z_0}{4} \left[ \partial_t p_0(t) - \partial_t p_0(t-T) \right]^2,
\end{equation}
where the left-moving current amplitude is a sum of the current emitted by the transmon and the delayed current arriving from the mirror.

In Fig.~\ref{fig:Energies}, we plot these energies for $Z_0/Z_J=1000$ for the case of resonance between transmon and the first cavity resonance $\omega_J=\omega_c$. The system energies indeed perform damped Rabi-oscillations with the frequency $\Omega=2/\sqrt{T C_J Z_0}$ and half the off-resonance damping rate $\gamma_J^m/2=1/2C_J Z0$, as indicated by the yellow line given by the expression $e^{-\gamma^m_J t/2} \cos^2 (\Omega t/2)$, approximating the full numerical solution of the differential equations very well. We note that Laplace transforming the equations of motion \eqref{eq:EoMphiJa}-\eqref{eq:V} and calculating the residues of the system variables, gives similar expressions for the Rabi frequency and damping rate as the analysis of the resonances in the scattering amplitudes. More details and a comparison of the approximation to the numerical results can be seen in the supplemental materials.

\section{Effective quantum model: atom in a multimode cavity Hamiltonian}
We now go on to demonstrate that in the high impedance regime the response function $f(\omega)$ of the field trapped between the transmon and the mirror Eq.~\eqref{f_omega}, reproduces the dynamics of an effective Hamiltonian of a single transition atom in a multimode cavity. This Hamiltonian is of the form:
\begin{align}
\label{Hn}
H = \omega_J\left(a^\+a +{\textstyle\frac{1}{2}}\right) &+ \sum_{n=1}^\infty n \omega_c\left(c_n^\+c_n + {\textstyle\frac{1}{2}}\right)  \nonumber \\ &+ \frac{\Omega}{2} \sum_{n=1}^\infty  \left(a^\+ + a\right)\left(c_n^\+ + c_n\right),
\end{align}
where the operators $a,a^\+$ are annihilation and creation bosonic operators ($[a,a^\+] = 1$) associated with excitation in the transmon qubit, while $c_n,c_n^\+$ annihilate/create photons in the cavity modes. When weakly excited, the choice for bosonic excitations of the transmon is justified, while the orthogonality relations between the cavity modes is ensured by the high finesse of the latter, so that we have $[c_n,c_m^\+] = \delta_{nm}$. Details of the diagonalisation of the Hamiltonian are shown in the supplemental materials. The response function $|f(\omega)|$ and eigenfrequencies of Hamiltonian \eqref{Hn} are shown in Fig.~\ref{eigenvalues}. The eigenfrequencies are shown to match the peaks of $|f(\omega)|$ for all cavity modes. Noticeably, a dip in the response function corresponding to the dark state is found for $\omega = \omega_0$.
\begin{figure}
	\centering
	\begin{overpic}[width=1\linewidth]{f_omega_v62.pdf}
		\put(98,21.53){3}
		\put(98,34.5){2}
		\put(98,47){1}
		\put(97,59.4){0}
		\put(97,72){1}
		\put(97,84.5){2}
	\end{overpic}
	\caption{Response function (log color scale) versus normalized frequency $\omega/\omega_J$ and linearly varied time delay  $\omega_JT$. The eigenvalues $\Omega_\alpha$ are computed for $N = 8$ cavity modes and are shown as superimposed white dashed lines. The parameters are chosen to be: $Z_0/Z_J = 100$ and $C_c/C_J = 0.3$.}
	\label{eigenvalues}
\end{figure}
\section{Discussion and Outlook}
We have made a first theoretical investigation of the properties of a transmon capacitively coupled to a high impedance transmission line, a system which is currently becoming experimentally accessible. By linearizing the Josephson junction, we could describe the low excitation dynamics, including spontaneous emission. We find qualitatively different behaviour, compared to the low impedance regime. In particular, the atom now forms its own cavity, and we can observe a vacuum Rabi splitting, giving rise to Rabi-oscillations in the spontaneous emission. The system is well described by a Hamiltonian for an atom weakly coupled to a multimode cavity.  We hope that this analysis will inspire an experimental realization of this novel system.

\section{Acknowledgements}
The authors acknowledge funding from the Swedish Research Council (VR) through Grant No. 2016-06059. GJ also acknowledge funding from the Knut and Alice Wallenberg foundation (KAW) through the Wallenberg Centre for Quantum Technology (WACQT).

\bibliographystyle{apsrev4-2}

\section*{Supporting information for Transmon in a semi-infinite high-impedance transmission line \textemdash\ appearance of cavity modes and Rabi oscillations}
\begin{center}
E. Wiegand$^1$, B. Rousseaux$^2$ and G. Johansson$^1$\\
\emph{$^1$Department of Microtechnology and Nanoscience - MC2, Chalmers University of Technology, 412 96 G\"oteborg, Sweden}\\
\emph{$^2$Laboratoire de Physique de l'\'Ecole Normale Sup\'erieure, ENS, Universit\'e PSL, CNRS, Sorbonne Universit\'e, Universit\'e de Paris, F-75005 Paris, France}
\end{center}

\author{Emely Wiegand}

\affiliation{Department of Microtechnology and Nanoscience - MC2, Chalmers University of Technology, 412 96 G\"oteborg, Sweden}

\author{Benjamin Rousseaux}
\affiliation{Laboratoire de Physique de l'\'Ecole Normale Sup\'erieure, ENS, Universit\'e PSL, CNRS, Sorbonne Universit\'e, Universit\'e de Paris, F-75005 Paris, France}

\author{G\"oran Johansson}
\affiliation{Department of Microtechnology and Nanoscience - MC2, Chalmers University of Technology, 412 96 G\"oteborg, Sweden}

\maketitle
\vspace{.5cm}

\section{Rabi oscillations of an atom coupled to a high-impedance semi-infinite TL - additional parameter values}
 We can Laplace transform the equations of motion of the transmon shown in the main article and extract the exact poles numerically. With the following formulas we can calculate the inverse Laplace transform of the system variables and energy:
\begin{align}
\text{Res}_{1,2}^{\pm} p_J(s) = \lim_{s \rightarrow s_{1,2}^{\pm}} p_J(s) (s - s_{1,2}^{\pm}),
\end{align}
\begin{align}
    p_J(t) = 	\int p_J(s) e^{-s t} d s =  \left.2 \pi i \sum_{k} \operatorname{Res} p_J(s) e^{-s t} \right|_{s=s_{k}}
    \label{eq:Laplace}
\end{align}
where $k = s_{1,2}^{\pm}$ are the poles of $p_J(s)$. Similarly, we calculate $\phi_J(t)$. We show the results as an addition to Fig.~4 in the main article. Here we provide more figures for different system parameters.
In all panels of Fig~\ref{fig:Rabi}, the impedance is chosen to be $Z_0/Z_J = 1000$. In Fig.~\ref{fig:Rabi} a), the ratio of the coupling capacitance and the Josephson capacitance is fairly small $C_c/C_J = 0.02$ and the coupling to the TL is weak. The cavity frequency equals the resonance frequency of the transmon $\omega_C= \omega_J$. The decay is weak and the Rabi-oscillations are clearly visible.  The parameters in Fig.~\ref{fig:Rabi} b) are similar to a), but now $C_c/(C_J+C_c) = 0.02$ and most importantly
the cavity frequency equals the transition frequency of the qubit for low impedance, $\omega_C= \omega_0$, which is the condition for a dark state \cite{Wiegand2020}. The energy of the transmon decays until it reaches the dark state, with energy $\frac{E_{DS}}{E_0} =  \frac{1}{(1+ \frac{T }{2}\gamma_0)^2}$ where $\gamma_0 = \frac{Z_0 \omega_0^2}{2}\frac{C_c^2}{C_c + C_J}$ \cite{Wiegand2020}. In c), the coupling capacitance is much larger compared to a) and b), $C_c/(C_c + C_J) = 0.3$ and the transmon fulfills the dark state condition
$\omega_C = \omega_0$. Anyhow, the system does not converge into a dark state and the Rabi oscillations are very weak. The main difference here is that, since $C_c/(C_c + C_J)$ is rather large, it means that $\omega_J$ is not close to $\omega_C$ and the Rabi oscillations and coupling to the cavity are suppressed. In this parameter regime, the behaviour of the system seems to be independent of the position of the transmon with respect to the mirror. Similar to Fig.~\ref{fig:Rabi} c), in Fig.~\ref{fig:Rabi} d), the coupling capacitance is rather large too $C_c/C_J = 0.3$, but here the cavity frequency equals the frequency of the transmon $\omega_C = \omega_J$, which means that the Rabi condition is fulfilled. We see clear Rabi oscillations and in addition the decay is much slower than in $c)$.
\begin{figure}[t]
	\flushleft
	\begin{minipage}{0.49\linewidth}
		\begin{overpic}[width=1\linewidth]{omegaJomegaC.pdf}
			\put(-4,60){ $a)$}
			\put(98,60){ $b)$}
			\put(-4,-8){ $c)$}
			\put(98,-8){ $d)$}
		\end{overpic}
	\end{minipage}
	\begin{minipage}{0.49\linewidth}
		\includegraphics[width=1\linewidth]{omega0omegaC.pdf}
	\end{minipage}
\flushleft
\begin{minipage}{0.49\linewidth}
	\includegraphics[width=1\linewidth]{Decay0.pdf}
\end{minipage}
\begin{minipage}{0.49\linewidth}
	\includegraphics[width=1\linewidth]{DecayJ.pdf}
\end{minipage}
\caption{Energy of the transmon (blue), energy of the field trapped between the transmon and the mirror (magenta), outging field to the left side of the transmon (orange), total energy of the system (red) and semi-numerically calculated energy of the transmon using Eq.~\eqref{eq:Laplace} as a function of period of the Rabi oscillations. The parameters for the panels are the following:
	a) $Z_0/Z_J = 1000$, $C_c/C_J = 0.02$ and $\omega_J = \omega_C$. Since we chose $C_c/C_J = 0.02$, the decay is very slow.
b) $Z_0/Z_J = 1000$, $C_c/(C_J+C_c) = 0.02$ and $\omega_0 = \omega_C$. Here, the parameters are almost the same as in a), but now the dark state condition is fulfiled, $\omega_J = \omega_0$. The system converges into a dark state with a finite excitation probability of both the transmon and the field between the transmon and the mirror. The green line indicates the analytical value of the dark state energy $\frac{E_{DS}}{E_0} =  \frac{1}{(1+ \frac{T }{2}\gamma_0)^2}$ with $\gamma_0 = \frac{Z_0 \omega_0^2}{2}\frac{C_c^2}{C_c + C_J}$.
c) $Z_0/Z_J = 1000$, $C_c/(C_J+C_c) = 0.3$ and $\omega_0 = \omega_C$. Here despite $\omega_0 = \omega_C$ as in b), the system does not seem to converge into a dark state and the decay rate is given by $\gamma_J = 1/C_J Z_0 $. The difference to b) is that the ratio between $\frac{C_c}{C_J} \frac{Z_0}{Z_J}$ here is much bigger than in b) which also means that $\omega_J$ is not close to $\omega_C$ and the Rabi oscillations are barely visible..
d) $Z_0/Z_J = 1000$, $C_c/C_J = 0.3$ and $\omega_J = \omega_C$. Here, as in c)  $\frac{C_c}{C_J} \frac{Z_0}{Z_J} \gg 1$  but the "Rabi condition" $\omega_J = \omega_C$ is still fulfilled. We see clear Rabi oscillations and the decay is slower compared to c) and the decay rate is given by $\gamma_J^m/2 = 1/2C_J Z_0 $.}
\label{fig:Rabi}
\end{figure}

As mentioned in the main article, we find an analytical expression for the oscillation frequency in the high impedance regime by analyzing the Laplace transform of the equations of motion. We find the Rabi frequency to be $\Omega = \frac{  2}{ \sqrt{T C_J Z_0} }$. In Fig.~\ref{fig:OmegaB}, we demonstrate the deviations of the approximation from numerically calculated values. We find that the higher the ratio $\frac{C_c}{C_J} \frac{Z_0}{Z_J}$, the closer the approximation resembles the numerical solution.
\begin{figure}[t]
	\centering
	\begin{overpic}[width=0.7\linewidth]{OmegaB2.pdf}
		\put(68,58){$\Omega^N (\frac{C_c}{C_J} = 0.01)$}
		\put(68,52){$\Omega^N (\frac{C_c}{C_J} = 0.05)$}
		\put(68,45){$\Omega^N (\frac{C_c}{C_J} = 0.3)$}
		\put(68,39){$\Omega^A$}
	\end{overpic}
	\caption{Frequency of the beating $\Omega$ as a function of $Z_0/Z_J$ calculated analytically (yellow, solid) and numerically (dotted) for different values of $C_c/C_J$, $C_c/C_J=0.01$ (purple), $C_c/C_J=0.05$ (red), $C_c/C_J=0.3$ (orange). We see that the higher $\frac{C_c}{C_J} \frac{Z_0}{Z_J}$ becomes, the closer the approximated analytical frequency $\Omega^A$ and the numerically calculated frequency $\Omega^N$ become.
	}
	\label{fig:OmegaB}
\end{figure}
\section{Analysis of the response functions}
To analyze the solution of the Fourier transformation of the equations of motion in the main article in a convenient manner, we introduce the following functions
\begin{align}
	N (\omega) &= R_0 (\omega) - i R_J (\omega), \\
	R_0 (\omega) &= \left( 1 - \frac{\omega^2}{\omega_0^2} \right) \\
	R_J (\omega) &= \frac{C_c  Z_0}{2} \omega \left( 1 - \frac{\omega^2}{\omega_J^2} \right).
\end{align}
With these definitions we are able to write the transmission and the reflection amplitudes of the qubit in the open transmission line as
\begin{align}
	r (\omega) = \frac{i R_J (\omega)}{ N (\omega)}, \quad
	t(\omega) = \frac{R_0 (\omega)}{N (\omega)}, \quad r(\omega) + 1 = t (\omega).
\end{align}
We note that the high impedance regime corresponds to that $|R_J(\omega)| \gg |R_0(\omega)|$ away from resonances, i.e. $|C_c Z_0 \omega /2 | \gg 1$, while the opposite ($|R_J(\omega)| \ll |R_0(\omega)|$) is true in the low impedance regime $|C_c Z_0 \omega /2 | \ll 1$.
\subsection{Damping rate for the open TL}
We analyse the scattering solution of the qubit excitation $\phi_J$ to find the damping rate for the transmon in an open TL. The solution in frequency space reads
\begin{align}
	\phi_J(\omega) = - \frac{C_c L_J \omega}{R_0 (\omega) - R_J (\omega)}.
\end{align}
In the high impedance regime, we perform an expansion around the bare qubit frequency $\omega = \omega_J + \delta \omega$ and find to first order in $\delta \omega$
\begin{align}
	\phi_J (\omega_J+\delta \omega) \approx \frac{i}{\omega_J}      \frac{1} {1 - i 2 \delta \omega/\gamma_J},
\end{align}
where $\gamma_J = 2/Z_0 C_J$ is the energy damping rate for spontaneous emission.

In the low impedance regime, the qubit resonance is shifted to $\omega_0$ and we instead expand $\omega = \omega_0 + \delta \omega$ to find
\begin{align}
	\phi_J (\omega_0+\delta \omega)  = \frac{1}{C_c Z_0 \omega_0^2} \frac{1}{1 - i 2 \delta \omega/ \gamma_0},
\end{align}
where $\gamma_0 = \frac{Z_0 \omega_0^2}{2} \frac{C_c^2}{C_c + C_J} $ is the low-impedance damping rate \cite{Wiegand2020}.
\subsection{Damping rates and Lamb shifts with a mirror}

With a mirror, the solution for the Josephson flux $\phi_J$
can be written as
\begin{align}
		\phi_J (\omega)  &=  \frac{  -i C_c L_J \omega \left(  1 - e^{i \omega T} \right) }{ R_0 (\omega) - i R_J (\omega)  \left(  1 - e^{ i \omega T} \right)  },
\end{align}
In the low impedance regime, we find that the qubit resonance is Lamb-shifted to $\tilde \omega_0 = \omega_0 + \gamma_0 \sin(\omega_0 T)/2$ and we can expand $\phi_J(\omega)$ around the resonance and find
\begin{align}
	\phi_J (\tilde{\omega}_0 +\delta \omega) \approx -i \frac{C_c}{C_J+C_c} \frac{e^{i\omega_0 T/2}}{\gamma_0 \sin \left( \omega_0 T/2\right)} \frac{1}{1 - i 2\delta \omega / \gamma_0^m},
\end{align}
with the Purcell-modified damping rate $\gamma_0^m = 2 \gamma_0 \sin^2 \left( \omega_0 T/2\right)$.

In the high impedance regime, the resonance frequency is Lamb shifted to $\tilde{\omega}_J=\omega_J+\gamma_J \cot{\left(\omega_J T /2\right)}/4$. Expanding $\phi_J(\omega)$ around this frequency we find
\begin{align}
	\phi_J (\tilde{\omega}_J + \delta \omega) = - \frac{2 i}{ \omega_J} \frac{1}{1 + i 2 \delta \omega /\gamma_J^m},
\end{align}
where $\gamma_J^m = \gamma_J/2=1/C_J Z_0$ is the damping rate of an atom in front of a mirror in the high $Z_0$ regime. Here, we note that the expression for the Lamb-shift diverges when $\sin{\left(\omega_J T /2\right)}=0$, i.e. when $\omega_J$ is close to a cavity resonance $\omega_C^n = 2 \pi n / T$. This is when the single pole approximation is no longer valid and we find the vacuum Rabi splitting. Away from the Rabi condition, we also note that the damping rate is independent of the distance to the mirror, i.e. we see no Purcell effect.
Away from the Rabi condition, we can also analyze the response function
\begin{align}
		f(\omega) &= \frac{1}{1 - i \frac{R_J (\omega)}{R_0 (\omega) }\left( 1 - e^{i \omega T}  \right)},
\end{align}
to extract the cavity modes.
Due to the finite transmission through the transmon, they are slightly shifted from the perfect mirror frequencies to $\tilde{\omega}_c^n=\omega_c^n + R_0(\omega_c^n)/T R_J(\omega_c^n)$. Close to the resonances we can expand $\omega=\tilde{\omega}_c^n+\delta \omega$ and find
\begin{align}
	f (\tilde{\omega}_c^n+\delta \omega) \approx \frac{2}{t (\omega_C^n)} \frac{1}{ 1 - i2\delta \omega / \gamma_c^n},
\end{align}
where $\gamma_c^n \approx |t(\tilde \omega_C^n)|^2/T$ is the energy damping rate of cavity mode $n$.

\section{Hopfield diagonalization of the atom in a mulimode cavity Hamiltonian}
To diagonalize the effective quantum optical Hamiltonian (Eq.~(12) in the main article), we first introduce the new polaritonic operators:
\begin{align}
	\Pi_\alpha = x^\alpha a + y^\alpha a^\+ +\sum_n\left(m_n^\alpha c_n + h_n^\alpha c_n^\+\right),
\end{align}
$x^\alpha,y^\alpha,m_n^\alpha,h_n^\alpha$ being the Hopfield coefficients associated to each bosonic operator $a,a^\+,c_n,c_n^\+$. To ensure the bosonicity of the polaritonic operators ($[\Pi_\alpha,\Pi_\beta^\+] = \delta_{\alpha\beta}$), the coefficients should satisfy the relation:
\begin{align}
	|x^\alpha|^2 - |y^\alpha|^2 + \sum_n\left(|m_n^\alpha|^2 - |h_n^\alpha|^2\right) = 1.
\end{align}
The new operators should satisfy the eigenvalue problem:
\begin{align}
	\left[\Pi_\alpha,H\right] = \Omega_\alpha\Pi_\alpha,
\end{align}
where $\Omega_\alpha$ are the eigenfrequencies labeled with a new index $\alpha$. Typically, if we couple the atom with $N$ cavity modes, then $\alpha$ runs from 1 to $N+1$. Expanding the commutator in the previous equation, it is possible to write the eigenvalue problem in a matrix form:
\begin{align}
	\left(
	\begin{array}{ccccccc}
	\omega_J & 0 & g_1 & -g_1 & \hdots & g_N & -g_N\\
	0 & -\omega_J & g_1 & -g_1 & \hdots & g_N & -g_N\\
	g_1 & -g_1 & \omega_1 & 0 & \hdots & 0 & 0\\
	g_1 & -g_1 & 0 & -\omega_1 &  & 0 & 0\\
	\vdots & \vdots & \vdots &  & \ddots & & \vdots \\
	g_N & -g_N & 0 & 0 &  & \omega_N & 0\\
	g_N & -g_N & 0 & 0 & \hdots & 0 & -\omega_N
	\end{array}\right)\left(\begin{array}{c}x^\alpha\\ y^\alpha\\m_1^\alpha\\h_1^\alpha\\\vdots\\m_N^\alpha\\h_N^\alpha\end{array}\right) = \Omega_\alpha\left(\begin{array}{c}x^\alpha\\ y^\alpha\\m_1^\alpha\\h_1^\alpha\\\vdots\\m_N^\alpha\\h_N^\alpha\end{array}\right).
\end{align}
This eigenvalue problem can be solved analytically for $N=1$ but in general, one has to diagonalize it numerically. A comparison between the full response function and the eigenvalues are shown in Fig.~5 in the main article.


\section{Introduction}

X-ray polarimeters are still at an early stage of development
where measurements are signal limited across the entire
band of sensitivity.
Current and planned instruments
generally provide sinusoidal signals
with fractional half-amplitudes of $p\mu$ for 100$p$\%
(linear) polarized input, where $\mu$
is the modulation factor.
Often, $p$ and $\mu$ are significantly
less than unity, providing weak signals.
The minimum detectable polarization (MDP) is frequently
used to indicate how low an instrument's sensitivity might be;
when the instrument gives $p > $ MDP,
the signal should be significant at high confidence.
Most calculations of MDP select a single value for $\mu$
in order to determine the MDP even
though $\mu$ may vary significantly with
energy within the instrument's bandpass.
Here, I examine how one may account for such
variation in an optimal fashion.

X-ray polarimeters that use the photoelectric effect
such as the Imaging X-ray Polarimetry Explorer \citep[IXPE:][]{2016SPIE.9905E..17W,ixpe} or
practically any other method
including scattering \citep[e.g. X-Calibur:][]{2014JAI.....340008B,10.1117/1.JATIS.4.1.011004},
provide a signal that is probabilistically related to the
input polarization direction.
A histogram of event phase angles, $\psi$, then has a characteristic
instrument-dependent
modulation factor, $\mu$.
Ground-based calibration using sources of known polarization angle are used to
determine how $\mu$ depends on energy $E$.
Polarimeters are usually characterized by their MDP which
depends on $\mu$ and is also related to the exposure time, $T$, and
the count rates the instrument records for both the source $R_S$ and background $R_B$:

\begin{equation}
{\rm MDP}_{99} = \frac{4.292}{\mu R_S} [\frac{R_S + R_B}{T}]^{1/2}
\end{equation}

\noindent
and the level of confidence, which we will take to be 99\% \citep{weisskopf:77320E}.
When background is negligible, then

\begin{equation}
\label{eq:mdp}
{\rm MDP}_{99} = \frac{4.292}{\mu N^{1/2}}
\end{equation}

\noindent
where $N = T R_S$ is the expected number of counts in the observation.

In general, the modulation factor depends on energy.
It is defined as the ratio of the semi-amplitude of the event
histogram to the average over phase for fully polarized light.
Thus, MDP$_{99}$ can be readily computed for a small energy range but it is more difficult
to determine what value of $\mu$ to use when $\mu(E)$ is a strong function
of $E$.  One approach would be to use a count-weighted average of $\mu$:

\begin{equation}
\label{eq:mdp0}
\bar{\mu} \equiv \frac{\sum_j \mu(E_j) C_j}{N}
\end{equation}

\noindent
where $C_j$ are the counts observed in energy bin $j$ and $N = \sum_j C_j$
\cite{2012SPIE.8443E..4NE}.
Here, I compute optimal estimates of the polarization Stokes parameters
for cases where $\mu$ can vary significantly across the energy range of the detector
as long as the Stokes parameters are constant or vary slowly across the entire range.

For high signal observations of strongly polarized sources, the data can be divided into
energy bands over which the assumption of $\mu = $ constant is appropriate.
In this analysis, I will concentrate on the case when the signal is weak, so that one
requires a large bandpass in order to detect a polarization signal.
As circular polarization is not yet feasibly detectable, I will assume that Stokes $V$ is zero.
Furthermore, for simplicity, only the case where the background is negligible and
the exposure for each angular bin is the same are considered.
\citet{2015APh....68...45K} provided
an method to determine weights for the case where an instrument
must be rotated to sample the Stokes parameters fully and the sampling is not uniform.

\section{Weighting of Binned Data}

\subsection{Assumptions}

Let $i$ designate one of $n$ angular phase bins of width $\Delta \psi = 2\pi/n$ about
$\psi_i = i \Delta \psi$ ($i = 0..n-1$) and $j$ indicate the
energy bin of width $\Delta E$ about $E_j$ ($j = 1..J$).  The instrument provides counts $C_{ij}$, each
with uncertainty $\sigma_{ij}$.  I assume that the background is negligible for this analysis.
Also, for the sake of simplicity so that we may gain a little intuition about
the solution, I will assume, for now, that we have large signals in every bin so that
we may use $\chi^2$ statistics.  To make formulae a little simpler, I will set $\sigma_{ij}$ to
$\sigma$ for a few of the derivations but then generalize in \S~\ref{sec:binned}.

Since we are interested in how to combine data from different energy channels, I also assume
that the Stokes parameters $Q$ and $U$ are not functions of energy.  Otherwise, we would
not combine the channels at all.  I will not use any formalism involving
energy to channel redistribution matrix files (aka RMFs), but instead assume
that the energy response functions are narrow compared to the energy bands.  This last
point is not really restrictive because all we need in the end is a model that describes
the counts in each phase-averaged energy bin.

\subsection{The Model and Fit Statistic}

The expected number of counts in bin $i,j$ is

\begin{equation}
\lambda_{ij} = [ 1 +  \mu_j  (q \cos 2\psi_i + u \sin 2\psi_i) ] f_j A_j T \Delta E \Delta \psi
\end{equation}

\noindent
where $f_j$ is the source flux in suitable units (say ph/cm$^2$/s/keV per radian of rotation), $\mu_j$ is the modulation
factor at energy $E_j$, $A_j$ is the system effective area, $T$ is the exposure time, and $q$ and $u$
are the fractional Stokes parameters, given by
$Q = P \cos \phi_0 = q I$ and
$U = P \sin \phi_0 = u I$, respectively.  The polarization fraction
is then $p \equiv P/I = (q^2 + u^2)^{1/2}$ and the source phase angle is $\phi_0 = \tan^{-1} u/q$.
The electric vector position angle (EVPA) is $\varphi = \phi_0/2$.

We form the fit statistic from the data and the model as
\begin{equation}
\chi^2 = \sum_i \sum_j \frac{(C_{ij} - \lambda_{ij})^2}{\sigma_{ij}^2}   .
\label{eq:chi2}
\end{equation}

\noindent
To make the solution a little easier to read, I define $\alpha_j = f_j A_j T$,
$\beta_j = \alpha_j \mu_j$, $s_i = \sin 2\psi_i$, and $c_i = \cos 2\psi_i$.

At this point, I set $\sigma_{ij} = \sigma$ for all $i$ and $j$, which would apply
if the counts in each energy bin were the same and $p \ll 1$.
Eq.~\ref{eq:chi2} then becomes

\begin{equation}
\chi^2 = \frac{1}{\sigma^2} \sum_i \sum_j (C_{ij} - [ \alpha_j + q \beta_j c_i + u \beta_j s_i ] \Delta E \Delta \psi)^2   .
\label{eq:chi2b}
\end{equation}

\noindent
Because the phase angles are chosen uniformly,

\begin{eqnarray}
\sum_i s_i = \sum_i c_i = 0 \\
\label{eq:crosssums}
\sum_i s_i  c_i = 0 \\
\sum_i s_i^2 = \sum_i c_i^2 = n/2  ~~~.
\label{eq:sumsq}
\end{eqnarray}

\subsection{Solving for the Flux (Stokes I)}

Setting the derivative of $\chi^2$ with respect to $\alpha_j$ to zero yields

\begin{eqnarray}
\sum_i C_{ij} & = & n \Delta' \alpha_j + \Delta' \mu_j^2 \alpha_j \sum_i (q^2 c_i^2 + u^2 s_i^2 ) \\
	& = & n \Delta' \alpha_j + \frac{n}{2} \Delta' \mu_j^2 p^2 \alpha_j
	\label{eq:exactalpha}
\end{eqnarray}

where $\Delta' = \Delta \psi \Delta E$.
For this analysis, $\mu_j p \ll 1$, so we drop the 2nd term on the right of Eq.~\ref{eq:exactalpha}, obtaining

\begin{equation}
\label{eq:alphaj}
\hat{\alpha}_j = \frac{\sum_i C_{ij}}{2\pi \Delta E} = \frac{C_{.j}}{2\pi \Delta E}
\end{equation}

\noindent
where $C_{.j} \equiv \sum_i C_{ij}$ is the expected number of counts in channel $j$.  We also have

\begin{equation}
\label{eq:fj}
\hat{f}_j = \frac{C_{.j}}{2\pi T  A_j \Delta E}   ~~~.
\end{equation}

\noindent
Note that the units are appropriate; the factor of $1/2\pi$ comes in because the flux was
defined per radian of polarization phase space.

\subsection{Solving for the Polarization (Stokes Q and U)}

Similarly proceeding, one may derive the best estimates of $q$ and $u$:
\begin{eqnarray}
\hat{q} & = & \frac{\sum_i c_i \sum_j \beta_j C_{ij} }{\Delta' \sum_i c_i^2 \sum_j \beta_j^2} \\
\hat{u} & = & \frac{\sum_i s_i \sum_j \beta_j C_{ij} }{\Delta'\sum_i s_i^2 \sum_j \beta_j^2} ~~~~.
\end{eqnarray}

\noindent
Using the definition of $\beta_j$ and substituting Eq.~\ref{eq:alphaj}, we have
\begin{eqnarray}
\label{eq:q}
\hat{q} = \frac{n \sum_i c_i \sum_j \mu_j  C_{.j} C_{ij}}{\sum_i c_i^2 \sum_j \mu_j^2 C_{.j}^2} =
	\frac{2 \sum_j \mu_j C_{.j} \sum_i c_i C_{ij}}{\sum_j \mu_j^2  C_{.j}^2} = \sum_i c_i x_i \\
\hat{u} = \frac{n \sum_i s_i \sum_j \mu_j C_{.j} C_{ij} }{\sum_i s_i^2 \sum_j \mu_j^2  C_{.j}^2} =
	\frac{2 \sum_j \mu_j C_{.j} \sum_i s_i C_{ij}}{\sum_j \mu_j^2  C_{.j}^2} = \sum_i s_i x_i
\label{eq:u}
\end{eqnarray}

\noindent
where the definition of $C_{.j}$ and Eq.~\ref{eq:sumsq} are used and

\begin{equation}
x_i \equiv \frac{2 \sum_j \mu_j C_{.j} C_{ij}}{\sum_j \mu_j^2 C_{.j}^2} ~~~.
\end{equation}

\noindent
Also,

\begin{equation}
\tan(\hat{\phi_0}) = \frac{\sum_j \mu_j C_{.j} \sum_i s_i C_{ij}}{\sum_j \mu_j C_{.j} \sum_i c_i C_{ij}} = \frac{\sum_i s_i x_{i}}{\sum_i c_i x_{i}}
\end{equation}

\noindent
\begin{equation}
\hat{p} = [ ( \sum_i c_i x_i )^2 + ( \sum_i s_i x_i )^2 ]^{1/2} = | \digamma_2(x) |
\end{equation}

\noindent
where $\digamma_2(x)$ is the 2nd term of the Fourier transform of $x_i$ (i.e., the $2\theta$ modulation
portion of the phase binned counts weighted by a channel-dependent term).
Note that $p \mu_j C_{.j}$ is the maximum modulated signal in channel $j$,  and $\mu_j C_{.j}$ is the maximum modulated signal when $p = 1$.

The best estimates of $q$ and $u$ for a single
energy channel $j$ are obtained from Eqs~\ref{eq:q} and \ref{eq:u} by removing the summation over $j$, giving

\begin{eqnarray}
\label{eq:q1}
\hat{q}_j = \frac{2 \sum_i c_i C_{ij}}{\mu_j  C_{.j} } \\
\hat{u}_j = \frac{2 \sum_i s_i C_{ij}}{\mu_j  C_{.j} }
\label{eq:u1}
\end{eqnarray}

\noindent
We can check these formulae
by modeling with $C_{ij} = C_{.j}/n (1 + \eta \mu_j \cos 2\psi_i ) = C_{.j}/n (1 + \eta \mu_j c_i )$ as the response of the detector to
a signal that is $100\eta$\% polarized at an EVPA $= \varphi = \phi_0 = 0$.
Substituting the model into Eqs.~\ref{eq:q1} and \ref{eq:u1}  gives
$\hat{q}_j = \eta$ and $\hat{u}_j = 0$, as expected, using Eqs.~\ref{eq:crosssums} and \ref{eq:sumsq}.

\subsection{Bandpass Weighting}

Solving Eqs.~\ref{eq:q1} and \ref{eq:u1} for $\sum_i c_i C_{ij}$ and $\sum_i s_i C_{ij}$, respectively,
then Eqs.~\ref{eq:q} and \ref{eq:u} can be rewritten as

\begin{eqnarray}
\hat{q} = \sum_j W_j \hat{q}_j \\
\hat{u} = \sum_j W_j \hat{u}_j
\end{eqnarray}

\noindent
where

\begin{equation}
W_j = \frac{\mu_j^2 C_{.j}^2}{\sum_k \mu_k^2 C_{.k}^2}
\end{equation}

\noindent
Because $\sum_j W_j = 1$, the $W_j$ values can be identified as weights to apply to the individual
bandpass values in order to obtain optimal, global values, as long as the polarization angle
and fraction do not depend significantly on energy.  The relative weight of channel $j$ compared
to channel $k$ is just $(\mu_j/\mu_k)^2 (C_{.j}/C_{.k})^2$, the ratio of the variances of the total
counts in the two channels weighted
by the square of the relevant modulation factors.

\subsection{Stokes Parameter Uncertainties, Binned Case, Constant Uncertainties}

\label{sec:constuncert}
The diagonal terms of the Fisher matrix give the parameter uncertainties for $f_j$, $q$, and $u$, as the
off-diagonal derivatives of Eq.~\ref{eq:chi2b} are all zero or very small when $p \mu \ll 1$.
Because of the assumption that the bandpasses are arranged to have equal
variances, $\sigma^2$, $C_{.j} = n \sigma^2 = N/J$, giving

\begin{eqnarray}
\sigma_{f_j}^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial f_j^2}} =
	\frac{2}{(T A_j)^2 \frac{ \partial^2 \chi^2}{\partial^2 \alpha_j^2}} = \frac{\sigma^2}{(T A_j \Delta')^2 n}
	= \frac{C_{.j}}{(2 \pi T A_j \Delta E)^2}\\
\sigma_q^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial q^2}} =
	\frac{2}{2 (\Delta')^2 \sum_i c_i^2 \sum_j \alpha_j^2 \mu_j^2} =
	\frac{2 n \sigma^2}{ \sum_j\mu_j^2 C_{.j}^2} = \frac{2 J}{N \sum_j \mu_j^2}  \\
\sigma_u^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial u^2}} =
	\frac{2}{2 (\Delta')^2 \sum_i s_i^2 \sum_j \alpha_j^2 \mu_j^2} =
	\frac{2 n \sigma^2}{ \sum_j\mu_j^2 C_{.j}^2} = \frac{2 J}{N \sum_j \mu_j^2} ~~~~.
\end{eqnarray}

\noindent
In order to relate to the Stokes parameters' uncertainties to MDP$_{99}$, consider
the case of a single channel, where $J = 1$ and  $\mu_j = \mu$, so

\begin{equation}
\sigma_q^2 = \sigma_u^2 = \frac{2}{\mu^2 N}~~~~.
\end{equation}

\noindent
MDP$_{99}$ is then just scaled to the uncertainty in $q$ or $u$:

\begin{equation}
{\rm MDP}_{99} = \frac{4.292}{\mu N^{1/2}} = 3.035 \sigma_q = 3.035 \sigma_u ~~~~.
\end{equation}

\noindent
Due to the symmetry between $q$ and $u$ in the absence of a signal, we can determine the
MDP for uniformly weighted polarimeter counts when there are several bandpasses with
differing modulation channels:

\begin{equation}
{\rm MDP}_{99} = 3.035 \sigma_q = \frac{4.292 J}{N \sum_j \mu_j^2} = \frac{4.292}{(\bar{\mu^2} N )^{1/2}} ~~~~.
\label{eq:mdp2}
\end{equation}

\noindent
Eq.~\ref{eq:mdp2} indicates that the optimal modulation factor to use
for a multi-band data set with constant counts per energy bin
is the RMS of the set of $\mu_j$ values, not the mean.

\subsection{Using Data Uncertainties, Binned Case}
\label{sec:binned}

While commonly practiced, it is not always appropriate to set spectral bins based on
the observed counts in order to obtain a constant number of counts per spectral bin.
One reason to avoid such an approach is that the bins may end up too wide to satisfy
the requirement that the modulation factor be constant across the bin.
So, now we return to the general case where uncertainties depend on the energy bin $j$.
With statistical weighting by $1/\sigma_{ij}^2$, the weighted sums over $c_i$, $s_i$, and $c_i s_i$
are no longer identically equal to zero as the uncertainties track the counts per bin, which will correlate
with $\psi_i$ for $p \mu_j > 0$.  However, when considering the special case of weakly
modulated signals, there can be simple solutions.

The best values of $\alpha_j$, $q$, and $u$ are found by setting $\nabla \chi^2 = 0$.
Using Eq.~\ref{eq:chi2} for $\chi^2$ gives

\begin{eqnarray}
\Delta' \hat{\alpha}_j  \sum_i (1 + \hat{q} \mu_j c_i + \hat{u} \mu_j s_i)^2 /\sigma_{ij}^2
	& = &  \sum_i C_{ij} (1 + \hat{q} \mu_j c_i + \hat{u} \mu_j s_i) /\sigma_{ij}^2, ~~j=1..J\\
\sum_j \mu_j^2 \hat{\alpha}_j^2 \sum_i (\hat{q} c_i^2 + \hat{u} c_i s_i)/ \sigma_{ij}^2 & = &
	\sum_j \mu_j \hat{\alpha}_j \sum_i (C_{ij}/\Delta' -   \hat{\alpha}_j) c_i / \sigma_{ij}^2 \\
\sum_j \mu_j^2 \hat{\alpha}_j^2 \sum_i (\hat{q} c_i s_i + \hat{u} s_i^2)/ \sigma_{ij}^2 & = &
	\sum_j \mu_j \hat{\alpha}_j \sum_i (C_{ij}/\Delta' -   \hat{\alpha}_j) s_i / \sigma_{ij}^2
\end{eqnarray}

\noindent
These equations will not be tractable in general but we can derive approximate solutions
for $p \mu_j \ll 1$ and we may drop the $c_i s_i$ terms that are small
compared to the $c_i^2$ and $s_i^2$ terms:

\begin{equation}
\hat{\alpha}_j = \frac{1}{\Delta'} \frac{\sum_i C_{ij} /\sigma_{ij}^2}{\sum_i 1 /\sigma_{ij}^2} \equiv \frac{\bar{C}_j}{\Delta'}
\end{equation}

\begin{eqnarray}
\label{eq:q2}
\hat{q} = \frac{\sum_j \mu_j \bar{C}_j \sum_i c_i (C_{ij}-\bar{C}_j) / \sigma_{ij}^2}{\sum_j \mu_j^2 \bar{C}_j^2 \sum_i c_i^2 / \sigma_{ij}^2}  \\
\hat{u} = \frac{\sum_j \mu_j \bar{C}_j \sum_i s_i (C_{ij}-\bar{C}_j) / \sigma_{ij}^2}{\sum_j \mu_j^2 \bar{C}_j^2 \sum_i s_i^2 / \sigma_{ij}^2}
\label{eq:u2}
\end{eqnarray}

\noindent
where the quantity $\bar{C}_j$ is a variance-weighted average of the counts in channel $j$.
When $\sigma_{ij}^2 = C_{ij}$,
then $\bar{C}_j = n/(\sum_i C_{ij}^{-1} )$.
The single channel estimates of the $q$ and $u$ values are

\begin{eqnarray}
\label{eq:qj2}
\hat{q_j} & = &  \frac{\sum_i c_i (C_{ij}-\bar{C}_j) / \sigma_{ij}^2}{\mu_j \bar{C}_j \sum_i c_i^2 / \sigma_{ij}^2}  \\
\hat{u_j} & = & \frac{\sum_i s_i (C_{ij}-\bar{C}_j) / \sigma_{ij}^2}{\mu_j \bar{C}_j \sum_i s_i^2 / \sigma_{ij}^2}  ~~~~.
\label{eq:uj2}
\end{eqnarray}

Again, we can solve Eqs.~\ref{eq:qj2} and \ref{eq:uj2} for the terms that sum over $C_{ij}/\sigma_{ij}^2$  in
the numerators to
obtain estimates of $q$ and $u$ in terms of single band estimates

\begin{eqnarray}
\label{eq:q3}
\hat{q} & = & \sum_j {w'}_j \hat{q}_j  \\
\hat{u} & = & \sum_j {w''}_j \hat{u}_j
\label{eq:u3}
\end{eqnarray}

\noindent
where

\begin{eqnarray}
{w'}_j = \frac{\mu_j^2 \bar{C}_j^2 \sum_i c_i^2/\sigma_{ij}^2}{\sum_k \mu_k^2 \bar{C}_k^2 \sum_i c_i^2/\sigma_{ik}^2} \\
{w''}_j = \frac{\mu_j^2 \bar{C}_j^2 \sum_i s_i^2/\sigma_{ij}^2}{\sum_k \mu_k^2 \bar{C}_k^2 \sum_i s_i^2/\sigma_{ik}^2}
\end{eqnarray}

\noindent
and, as before, $\sum_j {w'}_j = \sum_j {w''}_j = 1$.
These are the weights to be applied to single band results to combine them optimally.

\subsection{Stokes Parameter Uncertainties, Binned Case}

Again, uncertainties in the derived values can be estimated under the assumption $p\mu \ll 1$
so that cross-terms in the Fisher matrix are negligible.  Then,

\begin{eqnarray}
\sigma_{f_j}^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial f_j^2}} = \frac{1}{(T A_j \Delta')^2 \sum_i 1/\sigma_{ij}^2} \\
\sigma_q^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial q^2}} = \frac{1}{\sum_j  \mu_j^2 \bar{C}_j^2 \sum_i c_i^2/ \sigma_{ij}^2} \\
\sigma_u^2 & = & \frac{2}{\frac{\partial^2 \chi^2}{\partial u^2}} = \frac{1}{\sum_j  \mu_j^2 \bar{C}_j^2 \sum_i s_i^2/ \sigma_{ij}^2}
\end{eqnarray}

\noindent
For weak modulations, the counts in channel $j$ are approximately independent of phase bin, so
we may approximate $\sigma_{ij}^2$ by $\sigma_j^2 / n$, where $\sigma_j^2 \approx C_{.j}$ is the variance
in total counts in channel $j$, $C_{.j}$.  Now $\bar{C}_j = C_{.j}/n$, giving

\begin{eqnarray}
\sigma_q^2 & = &  \frac{n}{\sum_j  \mu_j^2 C_{.j}^2/\sigma_j^2 \sum_i c_i^2} =
	\frac{2}{\sum_j  \mu_j^2 C_{.j}^2/\sigma_j^2} \approx \frac{2}{\sum_j  \mu_j^2 C_{.j}}\\
\sigma_u^2 & = & \frac{n}{\sum_j  \mu_j^2 C_{.j}^2/\sigma_j^2 \sum_i s_i^2} =
	\frac{2}{\sum_j  \mu_j^2 C_{.j}^2/\sigma_j^2} \approx \frac{2}{\sum_j  \mu_j^2 C_{.j}}
\end{eqnarray}

Following the computation of MDP$_{99}$ as in \S~\ref{sec:constuncert}, MDP$_{99}$ in this case is

\begin{eqnarray}
\label{eq:mdp3a}
{\rm MDP}_{99} & = & 3.035 \sigma_q = \frac{4.292}{( \sum_j  \mu_j^2 C_{.j}^2/\sigma_j^2)^{1/2}} \\
	& \approx & \frac{4.292}{(\sum_j  \mu_j^2 C_{.j})^{1/2}} ~~~~.
\label{eq:mdp3}
\end{eqnarray}

\noindent
Eq.~\ref{eq:mdp3} indicates that the modulation factor to use in Eq.~\ref{eq:mdp}
that is optimal for a multi-band data set is the count-weighted
RMS of the set of $\mu_j$ values, $\sum_j \mu_j^2 C_{.j} / N$, not the count-weighted mean
(Eq.~\ref{eq:mdp0}).  The choice to use Eq.~\ref{eq:mdp3a} or Eq.~\ref{eq:mdp3} depends on whether
the variances are dominated by the total counts in the energy band; if there is
significant non-Poisson variance, then Eq.~\ref{eq:mdp3a} is
will be somewhat more accurate.

\section{Weighting of Unbinned Data}

In general, an unbinned likelihood analysis is more appropriate
than the use of $\chi^2$ statistics for fitting X-ray event data
due to the Poisson nature of the counting statistics.
Furthermore, binning details often involves user discretion.
Though normally requiring numerical methods to determine
model parameters, the assumption of weak modulations makes it
possible to obtain good estimates.

\subsection{Likelihood Formulation}

Assume that there are $N$ events, with energies and instrument phases $(E_i, \psi_i)$.
At energy $E_i$, the modulation factor is $\mu_i$, the instrument effective area is $A_i$,
and the intrinsic source flux is $f_i = f(E_i)$ based on the spectral model of the source.
The event density in a differential energy-phase element $dE d\psi$ about $(E, \psi)$ is

\begin{equation}
\lambda(E, \psi) = [ 1 +  \mu_E  (q \cos 2\psi + u \sin 2\psi) ] f_E A_E T dE d\psi
\end{equation}

\noindent
and the log-likelihood for a Poisson probability distribution of events, $S = -2 \ln L$, is

\begin{eqnarray}
S & = & -2 \sum_i \ln \lambda(E_i, \psi_i) + 2 T \int f_E A_E dE  \int_0^{2\pi} [ 1 +  \mu(E)  (q \cos 2\psi + u \sin 2\psi) ]  d\psi \\
 & = & -2 \sum_i \ln f_i -2 \sum_i  \ln( 1 +  q \mu_i \cos 2\psi_i + u \mu_i \sin 2\psi_i) + 4 \pi T \int f_E A_E dE + {\rm constant}
\label{eq:unbin}
\end{eqnarray}

\noindent
Note that eq.~\ref{eq:unbin} has two terms
that depend only on the parameters of $f_E$ and the remaining
interesting term depends only on
the parameters of $q$ and $u$.  If there are no parameters in common, in the
case where $q$ and $u$ do not depend on the source flux and we are
not fitting (yet) for parameters of the spectral model, then the log-likelihood
for the polarization parameters is merely

\begin{equation}
S(q,u) =  -2 \sum_i  \ln ( 1 +  q \mu_i \cos 2\psi_i + u \mu_i \sin 2\psi_i)   .
\label{eq:like}
\end{equation}

\subsection{Estimating $q$ and $u$ for the Unbinned Case}

We can define $c_i = \mu_i \cos 2\psi_i$ and
$s_i = \mu_i \sin 2\psi_i$ (unlike the case for binned data).  Setting
$\partial S/\partial q = 0$ and $\partial S/\partial u = 0$
to find the best estimates of $q$ and $u$ gives

\begin{eqnarray}
\label{eq:gen1}
0 = \sum_i \frac{c_i}{1 + \hat{q}c_i + \hat{u}  s_i} = \sum_i w_i c_i\\
0 = \sum_i \frac{s_i}{1 + \hat{q}c_i + \hat{u}  s_i} = \sum_i w_i s_i
\label{eq:gen2}
\end{eqnarray}

\noindent
where $w_i \equiv (1 + \hat{q} c_i + \hat{u}  s_i)^{-1}$.  It can be shown
that $\sum_i w_i = N$.  These two equations apply under quite general
circumstances but require numerical solution.  However, for $\hat{q} \ll 1$ and $\hat{u}  \ll 1$,
then $w_i \approx 1 - \hat{q} c_i - \hat{u} s_i$, so eqs.~\ref{eq:gen1} and \ref{eq:gen2} become

\begin{eqnarray}
\sum c_i = \hat{u} \sum c_i s_i + \hat{q} \sum c_i^2\\
\sum s_i = \hat{u} \sum s_i^2 + \hat{q} \sum c_i s_i
\end{eqnarray}

\noindent
where all sums are implicitly over $i$.  If the source is polarized, then the phase angles are not random, so the sums
over the sine or sine-cosine terms
are not identically zero.
Fortunately, the pair of linear equations can be solved, giving

\begin{eqnarray}
\hat{q} =  \frac{\sum c_i \sum s_i^2 - \sum s_i \sum c_i s_i}{\sum s_i^2 \sum c_i^2 - (\sum c_i s_i)^2}\\
\hat{u} =  \frac{\sum s_i \sum c_i^2 - \sum c_i \sum c_i s_i}{\sum s_i^2 \sum c_i^2 - (\sum c_i s_i)^2}
\end{eqnarray}

\noindent
In keeping with the assumed approximations of small $q$ and $u$,
we expect the sine-cosine terms to be small compared to the sums over $s_i^2$
or $c_i^2$, giving

\begin{eqnarray}
\hat{q} \approx \frac{\sum c_i}{ \sum c_i^2}\\
\hat{u} \approx \frac{\sum s_i}{ \sum s_i^2}
\end{eqnarray}

It should be simple to compute these sums for all the events in an observation in order to
obtain a first estimate of the sizes of the fractional Stokes parameters, $q$ and $u$.
These estimates provide some accounting for the variation of the modulation factor by
using it as a weight on each phase term.
Furthermore, the solution indicates a simple way to estimate the dependences of $q$ and $u$
on $E$ by merely creating the sums over suitably small energy ranges so one may obtain
a non-parametric estimate of the $(q,u)$ energy dependence, assuming that the values are
small but measurable.

The RMF does not enter this analysis, making it relatively easy to examine a data set
in a preliminary fashion before attempting detailed fitting that might require an
iterative fitting method that accounts for energy redistribution.
Technically, assigning a value of $\mu$ to an event
requires knowledge of the event's true energy, which
is uncertain because of the probabilistic
mapping between the detector signal and energy.
However, as long as the modulation factor
varies slowly compared to the detector response function, this approach may still
be useful and give a valid impression as to how to devise a physical model.

\subsection{Uncertainties, Unbinned Case}

Finally, the uncertainties on $q$ and $u$ are estimated from the 2nd derivatives
of eq.~\ref{eq:like}:

\begin{eqnarray}
\frac{\partial^2 S}{\partial u^2} = 2 \sum w_i^2 s_i^2\\
\frac{\partial^2 S}{\partial q^2} = 2 \sum w_i^2 c_i^2\\
\frac{\partial^2 S}{\partial u \partial q} = 2 \sum w_i^2 c_i s_i
\end{eqnarray}

\noindent
and using $w_i \approx 1$ and $\sum w_i^2 c_i s_i \approx 0$,
one obtains simple estimates for the Stokes uncertainties:

\begin{eqnarray}
\sigma_q^2 \approx \frac{2}{\frac{\partial^2 S}{\partial q^2}} \approx \frac{1}{\sum c_i^2}\\
\sigma_u^2 \approx \frac{2}{\frac{\partial^2 S}{\partial u^2}} \approx \frac{1}{\sum s_i^2}
\end{eqnarray}

\noindent
Remember that the $c_i$ and $s_i$ values have the modulation factor included
in this situation, so the denominators
are not $N/2$.  To estimate the denominators, we note that
\begin{eqnarray}
\sum c_i^2 + \sum s_i^2 = \sum \mu_i^2 = N \bar{\mu^2}  ~~,
\end{eqnarray}
where $\bar{\mu^2}$ is the
count-weighted average of $\mu^2$.
Furthermore, $c_i$ and $s_i$ should be similarly distributed
when $p \ll 1$, giving
\begin{eqnarray}
\sum c_i^2 \approx \sum s_i^2 \approx \frac{N}{2} \bar{\mu^2}
\end{eqnarray}

\noindent
Finally, as in \S~\ref{sec:binned}, we make the connection to MDP$_{99}$:

\begin{equation}
\label{eq:mdp4}
{\rm MDP}_{99} = 3.035 \sigma_q  \approx \frac{3.035}{(\sum_i c_i^2)^{1/2}}
	 \approx \frac{4.292}{(\bar{\mu^2} N)^{1/2}} ~~~~.
\end{equation}

\noindent
As in the binned case,
the modulation factor to use in Eq.~\ref{eq:mdp}
that is optimal for unbinned data set is the count-weighted
RMS of the set of $\mu_i$ values, not the count-weighted mean
(Eq.~\ref{eq:mdp0}).

\section{Example Application and Summary}

Armed with an appropriate weighting scheme, it is useful to determine how much
the MDP may improve using it.
To determine the effect of optimal binning, the results of an observations
can be computed using assumed functions for
the spectral model, effective area, and
$\mu(E)$.  It is important to get the dependences of these functions
on energy approximately correct; the normalizations are not relevant when just
selecting a total number of events for each simulated observation.

Fig.~\ref{fig:functions} shows assumed models of the effective area,
normalized to 1 at the maximum value, and of $\mu(E)$, both
for an instrument like IXPE \cite{ixpe}.
The standard modulation curve for IXPE uses a moment method
with ellipticity exclusion \citep{2003SPIE.4843..383B};
the data are from A.\ Di Marco, {\it et al.}, in preparation
\citep[see also][]{2016SPIE.9905E..17W}.
The effective area is approximate
\citep[as in][]{2012SPIE.8443E..4NE} -- only the shape is important for
this analysis.
For a given
total number of counts, $N$, two versions of MDP$_{99}$ were computed
using the area and modulation curves for each of two assumed spectral models.
The spectral model is a power law with spectral index $\Gamma = 0$ or 2,
where $f_E \propto E^{-\Gamma}$.
For case 1, MDP$_{99}$ is determined using a count weighted modulation
factor given by Eq.~\ref{eq:mdp0}.  For case 2, MDP$_{99}$ is given by
Eq.~\ref{eq:mdp3}.  Thus:

\begin{equation}
{\rm MDP}_{99} =  \frac{4.292 [\int A(E) E^{-\Gamma} dE]^{1/2} }{ [ N \int \mu(E)^X A(E) E^{-\Gamma} dE]^{1/2} } ~~~~.
\label{eq:mdp99}
\end{equation}

\noindent
where $X$ = 1 for case 1 and 2 for case 2 and the integral in the numerator normalizes the
spectrum so that there are $N$ expected counts.
For this exercise, $N = 10^5$ after the ellipticity exclusion,
achievable with \mbox{\em IXPE\/}\ for
sources with 2-8 keV fluxes of $\sim 10^{-11}$ erg cm$^{-2}$ s$^{-1}$ in
about 10 days \citep{ixpe}.
For Case 1, MDP$_{99}$ came out to 3.61\% and 4.52\%
for $\Gamma = 0$ or $2$, respectively.  For the optimum weighting of Case 2, MDP$_{99}$ values were
3.41\% and 4.21\%.  Thus, the improvement to the MDPs were 6.0\% and 7.2\% for the two
different spectral slopes. It is important to note that the fractional improvement is independent of $N$.

One can obtain further improvement in the MDP by increasing the modulation factor.
\citet{PEIRSON2021164740} show that an approach to measuring \mbox{\em IXPE\/}\ tracks using
convolutional neural networks can increase $\mu(E)$, as shown in Fig.~\ref{fig:functions},
and thereby improve the subsequent
estimate of the MDP.
The actual improvement in MDP in this approach depends somewhat on a weighting parameter, $\lambda$
\citep[cf.][]{PEIRSON2021164740}.
Using Eq.~\ref{eq:mdp99} with the $\lambda = 2$ modulation
factor curve
from \citet{PEIRSON2021164740}  and
assuming $N = 10^5$ can yield further reductions of
MDP$_{99}$ by 15\% and 8\% for the two spectral shapes to 2.97\% and 3.89\%, respectively.
For $\lambda = 1$, the improvement is less than a few percent.
\citet{PEIRSON2021164740} point out that the track uncertainty
weight method has no ellipticity exclusion, thus starting with more
events than with the moments method,
but does reduce the effective number of events by 15-20\%.

Clearly, there is an advantage to optimal weighting using the count-weighted
RMS of the modulation factor.  The actual improvement will depend on the details
of the shape of the effective area and modulation factor curves as well as the actual shape
of the spectrum.  However, this weighting can be generally applied for the best results.

 \begin{figure}
    \centering
    \includegraphics[width=15cm]{area_modulation_v2.pdf}
 \caption{
Representative functions for the effective area, $A$, and modulation
  function, $\mu$ for an X-ray polarimeter.  The effective area curve is illustrative,
  for an instrument like IXPE \citep{2012SPIE.8443E..4NE,ixpe}.  The
  effective area is normalized arbitrarily to the maximum value.
  The solid line gives the standard $\mu(E)$ for IXPE
  (A.\ Di Marco, {\it et al.}, in prep.).  The dashed curve gives the $\mu(E)$ curve
  from a neural network method with weighting parameter $\lambda = 2$ \cite{PEIRSON2021164740}.
 }
\label{fig:functions}
\end{figure}

\acknowledgments
Funding for this work was provided in part by contract 80MSFC17C0012 from the
Marshall Space Flight Center (MSFC) to MIT in support of IXPE, a NASA Astrophysics Small Explorers mission.
Support for this work was also provided in part by the National Aeronautics and
Space Administration (NASA) through the Smithsonian Astrophysical Observatory (SAO)
contract SV3-73016 to MIT for support of the Chandra X-Ray Center (CXC),
which is operated by SAO for and on behalf of NASA under contract NAS8-03060.
The support of the IXPE instrument Italian team, composed of
INAF and INFN personnel, is gratefully acknowledged.
I thank the anonymous referee for insightful comments that resulted in an improved manuscript.

\bibliographystyle{aasjournal}
\input{weighting_data.bbl}

\end{document}