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

Effective field theory (EFT) is an established tool in theoretical nuclear 
physics to connect the description of nuclei in terms of hadronic degrees of 
freedom to the underlying physics of quantum chromodynamics (QCD).  In relying 
only on basic concepts---symmetries, the separation of scales, and a systematic 
ordering of contributions---it is both elegant and pragmatic at the same 
time.  While the most widely used framework used for the description of nuclear 
structure and reactions is based on the approximate (and spontaneously 
broken) chiral symmetry of the QCD, for momentum scales well below the 
pion mass $M_\pi$, pion-exchange contributions cannot be resolved explicitly.  
The resulting pionless EFT contains, besides electromagnetic forces, only 
contact interactions between nonrelativistic
fields~\cite{Bedaque:1997qi,vanKolck:1997ut,Kaplan:1998tg,Bedaque:1998mb,
Kaplan:1998we,Birse:1998dk,vanKolck:1998bw,Chen:1999tn,Bedaque:1999vb,
Gabbiani:1999yv}.  From the original chiral symmetry it only keeps the isospin 
subgroup as an approximate symmetry.

This, however, is not the primary driving feature of the theory.  Instead, the 
fact that the $NN$ S-wave scattering lengths---$5.4$ (${-}23.7$)~\ensuremath{\mathrm{fm}}\ in 
the \ensuremath{{}^3S_1}\xspace (\ensuremath{{}^1S_0}\xspace) channel---happen to be large compared to the typical 
scale $M_\pi^{-1} \sim 1.4~\ensuremath{\mathrm{fm}}$ implies that the few-nucleon sector is close to 
a universal regime where short-range details have little impact on low-energy 
parameters.  In the two-nucleon sector, this is manifest in the effective range 
expansion (ERE) working as well as it does~\cite{Bethe:1949yr}.  In the 
three-nucleon sector, the spin-doublet S-wave configuration is governed 
by a non-derivative three-body interaction at leading order
(\text{LO}\xspace)~\cite{Bedaque:1999ve,Hammer:2000nf,Hammer:2001gh,Bedaque:2002yg,
Afnan:2003bs,Griesshammer:2005ga}, and the theory, formally equivalent to 
one for three bosons with short-range interactions, describes the triton 
as an approximate Efimov state~\cite{Efimov:1970zz,Efimov:1981aa,Hammer:2010kp}. 
For a recent review of the pionless three-body calculations, see 
Ref.~\cite{Vanasse:2016jtc}.

An important feature of the theory is that it can be constructed in a way such 
that it is fully renormalized order by order, with all 
nonperturbative effects included at \text{LO}\xspace and the rest treated in perturbation 
theory~\cite{Hammer:2001gh,Vanasse:2013sda,Konig:2013cia,Vanasse:2014kxa,
Konig:2015aka}.  This allows for precise calculations of low-energy 
observables with controlled error estimates based directly on the EFT expansion 
and independent of working at a fixed or limited regulator scale.  Recent 
examples of high-order calculations can be found in 
Refs.~\cite{Margaryan:2015rzg,Vanasse:2015fph}).  Moreover, pionless EFT 
provides a controlled ``laboratory'' to study such perturbative schemes 
(expanding directly the amplitude instead of a potential), which have 
been argued to be important for chiral EFT as well (see for example 
Refs.~\cite{Long:2016ybt,Valderrama:2016koj} for recent reviews).
  
While electromagnetic effects are certainly important for a realistic 
description of nuclei, their inclusion into the EFT provides a challenge 
because the long-range nature of these forces---meaning that the Coulomb force 
becomes dominant at very low energies, precisely where the short-range 
expansion otherwise works best---is not easily accommodated in the power 
counting.  Early work on Coulomb effects in the 
$pp$~\cite{Kong:1998sx,Kong:1999sf,Kong:2000px,Butler:2001jj,Barford:2002je,
Ando:2007fh,Ando:2008va,Ando:2008jb} and $pd$~\cite{Rupak:2001ci} systems was 
followed up on in Refs.~\cite{Ando:2010wq,Konig:2011yq}, with the latter 
providing a first calculation of $pd$ doublet-channel scattering.

It has become clear that if perturbative renormalization is to be maintained, 
including electromagnetic effects is not as simple as adding a Coulomb 
potential to the short-range terms, as it is typically done in calculations 
based on effective pionless potentials~\cite{Kirscher:2009aj,Kirscher:2011uc,
Kirscher:2011zn,Lensky:2016djr}.  Based on studying the regulator (cutoff) 
dependence of the amplitude, it was realized that in the presence of 
nonperturbative Coulomb effects an isospin breaking three-nucleon force is 
required to ensure renormalization at next-to-leading order
(\text{NLO}\xspace)~\cite{Vanasse:2014kxa,Konig:2014ufa}.  While another 
calculation~\cite{Kirscher:2015zoa} does not see the need to include such a 
term, Ref.~\cite{Konig:2015aka} showed that in part this three-body force is 
related to a two-body divergence in the $pp$ sector, and that the theory can be 
rearranged in such a way that all Coulomb effects in the \ensuremath{{}^3\mathrm{He}}\xspace bound state 
are included in perturbation theory.  In this counting the \ensuremath{{}^1S_0}\xspace channel is 
taken in the unitarity limit (infinite scattering length) at \text{LO}\xspace and the effects 
of the finite physical scattering lengths---$a_t$ as given above as well as the 
Coulomb-modified $pp$ scattering $a_C \simeq 7.81~\ensuremath{\mathrm{fm}}$---are accounted for by 
parameters entering at \text{NLO}\xspace.  In particular, in the $pp$ channel this parameter 
absorbs the logarithmic divergence generated by one-photon exchange in the 
two-nucleon subsystem.  Demoting all other (small) isospin-breaking effects to 
next-to-next-to-leading order (\text{N$^2$LO}\xspace) or higher, Ref.~\cite{Konig:2015aka} was 
able to calculate the \Triton--\ensuremath{{}^3\mathrm{He}}\xspace binding-energy difference at \text{NLO}\xspace 
without a new three-nucleon force.  More generally, this scheme enhances the 
predictive power of the theory by making the \ensuremath{{}^1S_0}\xspace channel parameter free at 
\text{LO}\xspace.

More recent work~\cite{Konig:2016utl} goes further and considers a more 
radical expansion that takes the \ensuremath{{}^3S_1}\xspace scattering length to infinity as well 
in a ``full unitarity'' leading order, whereas Ref.~\cite{Vanasse:2016umz} 
explores similar ideas by expanding around a leading order that exhibits the 
$SU(4)_W$ spin-isospin symmetry~\cite{Wigner:1936dx} for finite scattering 
lengths.

In this work, the pionless unitarity-\text{LO}\xspace counting schemes developed in 
Refs.~\cite{Konig:2015aka,Konig:2016utl} are implemented up to \text{N$^2$LO}\xspace, including 
effects from two-photon exchange and other isospin-breaking corrections.  It 
thereby establishes the convergence of these expansions up to this order.  
Going beyond \text{NLO}\xspace (first-order perturbation theory), where all corrections 
enter linearly, is an important proof of principle.  Moreover, it is 
demonstrated explicitly that away from the zero-energy threshold it is possible 
to describe $pd$ scattering with fully perturbative Coulomb effects.  Finally, 
it establishes the presence of an \text{N$^2$LO}\xspace three-body force, to be fixed by a 
single $pd$ doublet datum, at \text{N$^2$LO}\xspace.

\medskip
After discussing the basic setup and contributions---most notably two-photon 
diagrams---in Sec.~\ref{sec:Formalism}, the perturbative calculation of 
three-body observables is described in Sec.~\ref{sec:ThreeBody}.  Results
are presented and discussed in Sec.~\ref{sec:Results}, followed by a conclusion 
and outlook in Sec.~\ref{sec:Conclusion}.  Some details left out from the main 
text are provided in an appendix.

\section{Formalism and building blocks}
\label{sec:Formalism}

This section collects the ingredients required to set up the \text{N$^2$LO}\xspace 
calculation, summarizing results from previous works as far as necessary to 
make the current description reasonably self-contained.

\subsection{Effective Lagrangian and dibaryon propagators}

Using the notation of Ref.~\cite{Konig:2015aka}, the effective Lagrangian is 
split into one-, two- and three-body terms according to
\begin{equation}
 \mathcal{L} = N^\dagger\left(\mathrm{i} D_0+\frac{\boldsymbol{D}^2}{2M_N}\right)N
 +\mathcal{L}_\mathrm{2d}+\mathcal{L}_\mathrm{2t}+\mathcal{L}_3 \,,
\label{eq:L-Nd}
\end{equation}
where $N$ is the nonrelativistic nucleon field (doublet in spin and isospin 
space), coupled to (Coulomb) photons with the covariant derivative $D_\mu = 
\partial_\mu + \mathrm{i} eA_\mu \hat{Q}_N$ (with charge operator 
$\hat{Q}_N=(1+\tau_3)/2$).  The $NN$ S-wave two-body interactions are 
conveniently expressed in terms of auxiliary dibaryon fields $d^i$ (\ensuremath{{}^3S_1}\xspace) 
and $t^A$ (\ensuremath{{}^1S_0}\xspace), where $i$ and $A$ are spin-1 and isospin-1 indices, 
respectively.  The corresponding Lagrangians are
\begin{equation}
 \mathcal{L}_\mathrm{2d} = -d^{i\dagger}\left[\sigma_d
 + c_d\left(\mathrm{i} D_0+\frac{\boldsymbol{D}^2}{4M_N}\right)\right]d^i
 + y_d\left[d^{i\dagger}\left(N^T P^i_d N\right)+\mathrm{h.c.}\right] \,,
\label{eq:L-3S1}
\end{equation}
with the repeated superscripts $i$ summed over, and
\begin{multline}
 \mathcal{L}_\mathrm{2t} = -t^{0\dagger}\left[\sigma_t
 + c_t\left(\mathrm{i} D_0+\frac{\boldsymbol{D}^2}{4M_N}\right)\right]t^0
 - t^{{-1}\dagger}\left[\sigma_{t,pp}
 + c_{t,pp}\left(\mathrm{i} D_0+\frac{\boldsymbol{D}^2}{4M_N}\right)\right]t^{-1}\\
 - t^{{+1}\dagger}\left[\sigma_{t,nn}
 + c_{t,nn}\left(\mathrm{i} D_0+\frac{\boldsymbol{D}^2}{4M_N}\right)\right]t^{+1}
 +y_t\left[t^{{\tilde A}\dagger}
 \left(N^T \tilde{P}^{\tilde A}_t N\right)+\mathrm{h.c.}\right] \,.
\label{eq:L-1S0}
\end{multline}
Note that the \ensuremath{{}^1S_0}\xspace part is separated into physical channels ($np$, $nn$, 
$pp$).  Further details, including the projection operators $P^i_d$ and 
$\tilde{P}^{\tilde{A}}_t$ can be found in Ref.~\cite{Konig:2015aka}.  Setting
$y_d^2 = y_t^2 = 4\pi/M_N$, the remaining low-energy constants 
$\sigma_{d/t(,\cdot\cdot)}$ and $c_{d/t(,\cdot\cdot)}$ correspond directly to 
the scattering lengths and effective ranges in the respective channel.  
Contributions from the $NN$ shape parameters first enter at 
\text{N$^3$LO}\xspace~\cite{Margaryan:2015rzg} in the standard pionless power counting so that 
they also do not have to be included in the calculations presented here.

In the well-known fashion, nucleon-bubble insertions into the bare dibaryon 
propagators are resummed to obtain the full leading-order expressions
\begin{equation}
 \mathrm{i}\Delta_{d/t(,\cdot\cdot)}^{(0)}(p_0,\mathbf{p})
 = \frac{-\mathrm{i}}{\sigma_{d/t(,\cdot\cdot)}^{(0)}
 + y_{d/t(,\cdot\cdot)}^2 I_0(p_0,\mathbf{p})} \,,
\label{eq:Delta-dt}
\end{equation}
where
\begin{multline}
 I_0(p_0,\mathbf{p}) = M_N\int^\Lambda\dq{q}
 \frac{1}{M_N p_0 - \mathbf{p}^2/4 - \mathbf{q}^2 + \mathrm{i}\varepsilon} \\
 = {-}\frac{M_N}{4\pi}
 \left(\frac{2\Lambda}{\pi} - \sqrt{\frac{\mathbf{p}^2}{4}-M_N p_0-\mathrm{i}\varepsilon}\right)
 +\mathcal{O}(1/\Lambda)
\label{eq:I0-cutoff}
\end{multline}
is the generic nucleon bubble integral (Green's function from zero to zero 
separation) calculated with a sharp momentum cutoff.

In principle, an S-D mixing operator in the \ensuremath{{}^3S_1}\xspace channel enters at \text{N$^2$LO}\xspace. 
However, at this order it does not contribute to the \Triton--\ensuremath{{}^3\mathrm{He}}\xspace binding 
energy splitting, and neither is it relevant for the doublet S-wave 
phase shift~\cite{Vanasse:2013sda}.  Hence, it is not necessary to consider 
this operator in the present work.  Moreover, terms stemming from relativistic 
corrections and transverse photons are suppressed by inverse powers of $M_N$ 
and do thus not contribute to the order considered in this paper.

\subsubsection{Power counting around the unitarity limit}

In the power counting developed in Ref.~\cite{Konig:2015aka}, which is applied 
here up to \text{N$^2$LO}\xspace, the usual pionless expansion in terms of $Q/\Lambda_{\slashed\pi}$,
where $Q$ denotes the typical momentum of the process under 
consideration and $\Lambda_{\slashed\pi} \sim M_\pi$ is the pionless breakdown scale,
is paired with an additional expansion in $\aleph_0/Q$, where
\begin{equation}
 \aleph_0 \sim \alphaM_N \sim 1/|a_t| \sim 1/|a_C|
\end{equation}
is a typical low-energy scale in the \ensuremath{{}^1S_0}\xspace channel.  This combines the 
inverse scattering lengths with the typical Coulomb scale $\alphaM_N$.  In 
particular, at leading order the \ensuremath{{}^1S_0}\xspace channel is considered in the 
unitarity limit (where the scattering lengths are infinite) and thus also 
isospin symmetric.  The expansion in $\aleph_0/Q$ then corresponds to including 
the effects of finite $1/a_t$ and $1/a_C$ in perturbation theory, where the 
latter is naturally paired with one-photon exchange and ensures consistent 
renormalization in the presence of Coulomb effects.  Range corrections reflect 
the $Q/\Lambda_{\slashed\pi}$ expansion, as in the standard pionless counting.  Details of 
how the expansion is implemented by fixing the dibaryon parameters are given in 
Sec.~\ref{sec:TwoBodyParamaters}.

This scheme is constructed for the regime $\aleph_0 \ll Q \ll \Lambda_{\slashed\pi}$, which 
naturally includes the \ensuremath{{}^3\mathrm{He}}\xspace bound states.  This was studied in 
Ref.~\cite{Konig:2015aka}, which also found that the $\aleph_0/Q$ expansion
works well in the $nd$ scattering system.  In $pd$ scattering, Coulomb effects 
are certainly nonperturbative at very small center-of-mass momenta ($k\to0$), 
but for larger $k$ (determining $Q$ in this case), the perturbative expansion 
should work as well.  This is demonstrated in the present work.

More generally, the scheme of Ref.~\cite{Konig:2015aka} counts isospin-breaking 
corrections, not only of electromagnetic origin but also those induced by the 
up-down quark-mass difference in QCD.  As required by renormalization, the 
effects that give rise $a_C \neq a_t$ to are accounted for at \text{NLO}\xspace together 
with electromagnetic contributions.  Different values have been determined for 
the $nn$ scattering length~\cite{GonzalezTrotter:2006wz,Huhn:2001yk}, but it is 
generally assumed large (and negative), such that $1/\ensuremath{a_{t,nn}} \approx 1/a_t$.  
Following the counting of Ref.~\cite{Konig:2015aka}, this difference is 
accounted for here by an \text{N$^2$LO}\xspace correction, and for definiteness the central 
value favored by the pionless analysis of Ref.~\cite{Kirscher:2011zn}, 
$\ensuremath{a_{t,nn}} \simeq {-}(22.9\pm 4.1)~\ensuremath{\mathrm{fm}}$, is adopted here.  The small isospin 
breaking in the effective ranges, $r_C - r_t \simeq 0.06~\ensuremath{\mathrm{fm}}$ is also included 
at \text{N$^2$LO}\xspace.\footnote{The individual values used here are 
$r_C=2.794~\ensuremath{\mathrm{fm}}$~\cite{Bergervoet:1988zz} and 
$r_t=2.73~\ensuremath{\mathrm{fm}}$~\cite{Preston:1975}.}

\medskip
Recently, Ref.~\cite{Konig:2016utl} suggested a more radical rearrangement 
of the power counting that includes the \ensuremath{{}^3S_1}\xspace inverse scattering length 
$1/a_d$ in $\aleph_0$ instead of counting it as $\mathcal{O}(Q)$ (which is done in the 
standard pionless counting).  Although this expansion only perturbatively 
moves the deuteron bound state---which at \text{LO}\xspace is located at zero energy in this 
expansion---to its physical position, Ref.~\cite{Konig:2016utl} found that 
it works well for three and four-nucleon bound states at \text{NLO}\xspace.  As part of 
this work, \ensuremath{{}^3\mathrm{He}}\xspace is considered up to second order in this ``full unitarity'' 
scheme.

\subsubsection{Three-nucleon forces}

It is well known that in pionless EFT a three-nucleon interaction is needed 
already at \text{LO}\xspace to ensure renormalization of the doublet-channel 
amplitude~\cite{Bedaque:1999ve}.  This piece can be written
as~\cite{Ando:2010wq,Griesshammer:2011md}
\begin{equation}
 \mathcal{L}_{3,\text{LO}\xspace} = \frac{h}{3}N^\dagger\left[y_d^2\,
 d^{i\dagger} d^j \sigma^i \sigma^j+y_t^2\,t^{A\dagger} t^B \tau^A\tau^B
 - y_dy_t\left(d^{i\dagger} t^A \sigma^i \tau^A + \mathrm{h.c.}\right) \right]N \,,
\label{eq:L-3}
\end{equation}
where the coupling $h$ is also split up in various orders, with each piece 
having a characteristic log-periodic dependence on the momentum cutoff 
$\Lambda$ that is used in the three-body integral equations:
\begin{equation}
 h = \frac{2H_{0,0}(\Lambda)}{\Lambda^2}
 + \frac{2H_{0,1}(\Lambda)}{\Lambda^2}
 + \frac{2H_{0,2}(\Lambda)}{\Lambda^2}
 + \cdots \,.
\end{equation}
In addition to this, Ref.~\cite{Ji:2012nj} firmly established (for the 
analogous three-boson system) that an additional energy-dependent three-body 
force enters at \text{N$^2$LO}\xspace.  Following Ref.~\cite{Vanasse:2013sda}, where this was 
worked out for the $nd$ doublet system, it is included here in the form
\begin{equation}
 h_2 = \frac{2H_{2,2}(\Lambda)}{\Lambda^4}(M_N E + \gamma_d^2) \,,
\end{equation}
which conveniently vanishes at the $Nd$ scattering threshold and thus 
simplifies the numerical parameter determination.  These standard three-body
force are not changed by expanding around the unitarity limit or the 
perturbative inclusion of Coulomb effects.

\medskip
While the \text{NLO}\xspace isospin-breaking three-nucleon force identified 
in~\cite{Vanasse:2014kxa} is not needed in the counting scheme employed here 
(due to the fully perturbative treatment of Coulomb effects and only using 
isospin-symmetric effective ranges in the \ensuremath{{}^1S_0}\xspace channel at that order), it 
turns out (see Sec.~\ref{sec:ThreeBodyEquations}) that such a term is 
eventually needed at \text{N$^2$LO}\xspace.  In Sec.~\ref{sec:Results} it is shown that a term 
of the form
\begin{equation}
 h_\alpha = \frac{2H_{0,2}^{(\alpha)}(\Lambda)}{\Lambda^2} \,,
\end{equation}
with an operator structure designed to act only in the $pd$ doublet channel but 
otherwise give the same factors as for $h$~\cite{Vanasse:2014kxa}, is 
sufficient to ensure renormalization at \text{N$^2$LO}\xspace.

\subsection{Second-order binding energy shifts}
\label{sec:SecondOrderEnergies}

The off-shell amplitude, determined by the Lippmann--Schwinger equation 
describing nucleon-deuteron scattering (see Sec.~\ref{sec:ThreeBodyEquations}), 
is the central object considered in this paper.  It describes both scattering 
(in the on-shell limit) as well as bound-state properties.  The binding energy 
of a given state can be obtained from it by expanding the expression that 
formally includes all orders, as discussed in Ref.~\cite{Ji:2012nj}.  In the 
limit where the energy $E$ approaches the bound-state pole, it can be written as
(with $k$ and $p$ denoting in- and outgoing momenta, respectively)
\begin{multline}
 \mathcal{T}^{(0)}(E;k,p) + \mathcal{T}^{(1)}(E;k,p) + \mathcal{T}^{(2)}(E;k,p) + \cdots \\
 = {-}\frac{\mathcal{Z}^{(0)}(k,p) + \mathcal{Z}^{(1)}(k,p) + \mathcal{Z}^{(2)}(k,p) + \cdots}
 {E + B_0 + B_1 + B_2 + \cdots}
 + \mathcal{R}^{(0)}(E;k,p) + \cdots \,.
\label{eq:T-pole-012}
\end{multline}
In this expression, $\mathcal{Z}^{(n)}$ and $\mathcal{R}^{(n)}$, respectively, denote 
residue and regular terms, and the superscripts label contributions from 
different orders (note that regular terms with $n>0$ have been omitted in 
Eq.~\eqref{eq:T-pole-012}).  Expanding this by formally factoring out a small 
parameter and matching orders naturally recovers
\begin{equation}
 \mathcal{T}^{(0)}(E;k,p) = {-}\frac{\mathcal{Z}^{(0)}(k,p)}{E + B_0} + \mathcal{R}^{(0)}(E;k,p)
 \mathtext{with} \mathcal{Z}^{(0)}(k,p) = \mathcal{B}^{(0)}(k)\,\mathcal{B}^{(0)}(p)
\label{eq:T0-B0}
\end{equation}
at leading order.  Using the vertex functions $\mathcal{B}^{(0)}$, which can be 
obtained directly by solving a homogeneous integral equation and properly 
normalizing the solutions~\cite{Konig:2011yq}, one can write down the various 
matrix elements that contribute to the \text{NLO}\xspace energy shift.  This corresponds 
directly to first-order perturbation theory with bound-state wavefunctions and 
is what was done in Refs.~\cite{Konig:2014ufa,Konig:2015aka}.
Ref.~\cite{Ji:2012nj} instead starts from the \text{NLO}\xspace matching condition,
\begin{multline}
 \mathcal{T}^{(1)}(E;k,p)
 = B_1 \frac{\mathcal{Z}^{(0)}(k,p)}{(E + B_0)^2} + \frac{\mathcal{Z}^{(1)}(k,p)}{E + B_0}
 + \mathcal{R}^{(1)}(E;k,p) \\
 \implies B_1 = \lim\nolimits_{E\to{-}B_0} \frac{(E+B_0)^2\,\mathcal{T}^{(1)}(E;k,p)}
 {\mathcal{B}^{(0)}(k)\,\mathcal{B}^{(0)}(p)} \,,
\label{eq:B1-lim}
\end{multline}
and simplifies this to recover again the explicit expression in terms of 
$\mathcal{B}^{(0)}$.  At \text{N$^2$LO}\xspace,
\begin{equation}
 B_2 = \lim\nolimits_{E\to{-}B_0}
 \frac{(E+B_0)^2\,\mathcal{T}^{(2)}(E;k,p) + B_1(E+B_0)\mathcal{T}^{(1)}(E;k,p)}
 {\mathcal{B}^{(0)}(k)\,\mathcal{B}^{(0)}(p)} \,,
\label{eq:B2-lim}
\end{equation}
which Ref.~\cite{Ji:2012nj} expands explicitly in terms of $\mathcal{B}^{(0)}$ and 
$\mathcal{R}^{(0)}(E;k,p)$.  Since the expressions become rather involved at this 
point (particularly so with the inclusion of Coulomb corrections for \ensuremath{{}^3\mathrm{He}}\xspace), 
it is more convenient to directly evaluate Eq.~\eqref{eq:B2-lim}.  Numerical 
solutions at any given order can be obtained very efficiently using the 
technique described in Ref.~\cite{Vanasse:2013sda}.  In 
Ref.~\cite{Konig:2013cia}, it was extended to include Coulomb contributions and 
applied to calculate the quartet-channel $pd$ scattering up to \text{N$^2$LO}\xspace.

\subsection{Two-body T-matrix at second order}
\label{sec:TwoBody-NNLO}

To demonstrate how this procedure works in general, it is instructive to 
take as an example the \ensuremath{{}^3S_1}\xspace nucleon-nucleon system expanded around 
the unitarity limit, which was considered recently in Ref.~\cite{Konig:2016utl}. 
The calculation here also serves to prove explicitly what was claimed for the 
\text{N$^2$LO}\xspace deuteron energy shift in that paper.

Switching temporarily to a pionless formalism without dibaryon fields, the 
momentum-independent contact interaction $C_0$ is expanded as
\begin{equation}
 C_0 = C_0^{(0)} + C_0^{(1)} +  C_0^{(2)} + \cdots \,.
\label{eq:C0-expansion}
\end{equation}
No explicit label is included here to keep the notation simple because this 
subsection considers only a single channel.

\begin{figure}[htbp]
\centering
\begin{minipage}{0.95\textwidth}
\includegraphics[height=3.3em]{NN-IntEq-0}\\[0.66em]
\includegraphics[height=3.3em]{NN-IntEq-1}\\[0.66em]
\includegraphics[height=3.3em]{NN-IntEq-2}
\end{minipage}
\caption{Diagrammatic $NN$ integral equations up to second order.  Solid 
lines represent nucleons and an open circles with number $n$ denotes the $n$-th 
order T-matrix part.  A dot with $i$ surrounding circles corrections to a 
$C_0^{(i)}$ vertex.}
\label{fig:NN-IntEq-012}
\end{figure}

Furthermore ignoring range corrections (\textit{i.e.}, momentum-dependent contact terms) 
for simplicity, the integral equations that define the T-matrix up to second 
order in the approach of Ref.~\cite{Vanasse:2013sda} are shown in 
Fig.~\ref{fig:NN-IntEq-012}.  Note that lower-order amplitudes appear as part of 
the \text{NLO}\xspace and \text{N$^2$LO}\xspace inhomogeneous terms, but that the integral-equation 
kernel---determined solely by $C_0^{(0)}$---is always the same.  Iterating 
these equations shows that they recover the standard distorted-wave amplitudes 
that can also be obtained by a direct calculation of the individual diagrams at 
each order.  The primary advantage of the integral-equation formalism is that 
it avoids the explicit calculation of a two-loop integral (involving a full 
off-shell \text{LO}\xspace T-matrix), which was the primary motivation for its introduction 
in Ref.~\cite{Vanasse:2013sda}.\footnote{In the two-body sector with separable 
regularization, even full off-shell contributions become trivial and can be 
calculated analytically.  In the three-body sector, however, using the 
integral-equation technique is a substantial simplification.}

All two-nucleon bubbles in Fig.~\ref{fig:NN-IntEq-012} correspond to the 
integral $I_0$ given in Eq.~\eqref{eq:I0-cutoff}.  With the sharp cutoff 
regulator (or more generally any separable one), the integral equations at each 
order can be solved algebraically.  At leading order (first line in 
Fig.~\ref{fig:NN-IntEq-012}), one simply recovers the well-known result
\begin{equation}
 \mathcal{T}^{(0)}(E) = {-}C_0^{(0)} + C_0^{(0)} I_0(E)\,\mathcal{T}^{(0)}(E)
 = \frac{{-}C_0^{(0)}}{1 - C_0^{(0)}I_0(E)} \,,
\label{eq:T-NN-0}
\end{equation}
corresponding directly to the propagator expression~\eqref{eq:Delta-dt}.  At 
\text{NLO}\xspace and \text{N$^2$LO}\xspace, the equations and solutions are
\begin{equation}
 \mathcal{T}^{(1)}(z) = {-}C_0^{(1)} + C_0^{(1)} I_0(z)\,\mathcal{T}^{(0)}(z)
 + C_0^{(0)} I_0(z)\,\mathcal{T}^{(1)}(z) \\
 = \frac{{-}C_0^{(1)}\left[1 + I_0(z)\,\mathcal{T}^{(0)}(z)\right]}
 {1 - C_0^{(0)}I_0(z)}
\label{eq:T-NN-1}
\end{equation}
and
\begin{multline}
 \mathcal{T}^{(2)}(z) = {-}C_0^{(2)} + C_0^{(2)} I_0(z)\,\mathcal{T}^{(0)}(z)
 + C_0^{(1)} I_0(z)\,\mathcal{T}^{(1)}(z)
 + C_0^{(0)} I_0(z)\,\mathcal{T}^{(2)}(z) \\
 = \frac{{-}C_0^{(2)}\left[1 + I_0(z)\,\mathcal{T}^{(0)}(z)\right]
 - C_0^{(1)}I_0(z)\,\mathcal{T}^{(1)}(z)}{1 - C_0^{(0)}I_0(z)} \,.
\label{eq:tau-2}
\end{multline}

The leading-order term $C_0^{(0)}$, chosen to get a deuteron 
state with binding momentum $\kappa = \kappa^{(0)} \to 0$, corresponds to 
approaching the unitarity limit from the side of positive scattering length 
(\textit{i.e.}, a bound state moves towards zero energy).  At \text{NLO}\xspace, $C_0^{(1)}$ is 
adjusted to get the first term in the effective range expansion given by the 
physical scattering length, $\kappa^{(1)} = 1/a_d$.  These conditions are 
satisfied by setting
\begin{equation}
 C_0^{(0)} = {-}\frac{2\pi^2}{M_N\Lambda}
 \left(1 - \frac{\pi\kappa^{(0)}}{2\Lambda}\right)^{-1} \mathtext{,}
 C_0^{(1)} = {-}\frac{M_N\kappa^{(1)}}{4\pi} \big(C_0^{(0)}\big)^2 \,.
\label{eq:C00-C01}
\end{equation}
At \text{N$^2$LO}\xspace, they are maintained with
\begin{equation}
 C_0^{(2)}
 = \left(\frac{M_N\kappa^{(1)}}{4\pi}\right)^2 \big(C_0^{(0)}\big)^3 \,.
\label{eq:C02}
\end{equation}
Using a different regularization scheme (\textit{e.g.}, dimensional regularization with 
power divergence subtraction~\cite{Kaplan:1998tg} or the Gaussian regulator 
used for the four-body calculations in Ref.~\cite{Konig:2016utl}) only changes 
the details of $C_0^{(0)}$ but not the general features of the calculation.  
By construction, this gives
\begin{equation}
 \mathcal{T}^{(0)}(E)
 = \frac{4\pi}{M_N\left(\sqrt{\mathstrut{-}M_N E}-\kappa^{(0)}\right)}
 = \frac{8\pi\kappa^{(0)}}{M_N^2\left(E + \frac{(\kappa^{(0)})^2}{M_N}\right)}
 + \text{reg.\ terms}
\label{eq:T-0}
\end{equation}
and
\begin{equation}
 \mathcal{T}^{(1)}(E) = \frac{4\pi\kappa^{(1)}}
 {M_N\left(\sqrt{\mathstrut{-}M_N E}-\kappa^{(0)}\right)^{2}}
 \mathtext{,}
 \mathcal{T}^{(2)}(E) = \frac{4\pi\left(\kappa^{(1)}\right)^2}
 {M_N\left(\sqrt{\mathstrut{-}M_N E}-\kappa^{(0)}\right)^{3}} \,.
\label{eq:T-12}
\end{equation}
The \text{LO}\xspace binding energy and the shifts up to second order are found to be
\begin{equation}
 B_0 = \frac{(\kappa^{(0)})^2}{M_N}
 \mathtext{,}
 B_1 = \frac{2\kappa^{(0)}\kappa^{(1)}}{M_N}
 \mathtext{,}
 B_2 = \frac{(\kappa^{(1)})^2}{M_N}
\end{equation}
from Eqs.~\eqref{eq:T0-B0}, \eqref{eq:B1-lim}, and~\eqref{eq:B2-lim}.  In 
particular, in the unitarity limit $\kappa^{(0)}\to0$ the deuteron remains at 
zero energy at \text{NLO}\xspace, but it can be seen explicitly that it moves to its 
zero-range position, $B_2 \sim (\kappa^{(1)})^2 = 1/a_d^2$, at \text{N$^2$LO}\xspace,
exactly as claimed in Ref.~\cite{Konig:2016utl}.  While the results in 
Eqs.~\eqref{eq:T-0} and~\eqref{eq:T-12} match what one would naïvely expect 
from a direct expansion of the renormalized amplitude (\textit{i.e.}, written in terms of 
$\kappa = \kappa^{(0)} + \kappa^{(1)}$), it is reassuring to see the 
binding-energy shifts come out as desired even though in the limit 
$\kappa^{(0)}\to0$ the deuteron disappears as a state with normalizable 
wavefunction.

\subsection{Two-photon contributions}

As mentioned in the introduction, one motivation behind constructing the 
\ensuremath{{}^1S_0}\xspace unitarity-limit scheme of Ref.~\cite{Konig:2015aka} was that it makes 
it very convenient to include perturbative Coulomb corrections in a way that 
ensures proper renormalization.  Specifically, the log-divergent diagram, where 
a Coulomb photon is exchanged inside a $pp$ bubble, is included at \text{NLO}\xspace along 
with $\sigma_{t,pp}^{(1)}$, which absorbs the divergence when it is adjusted to 
give the physical (Coulomb-modified) $pp$ scattering length.  In this scheme, 
two-photon contributions (see Figs.~\ref{fig:TwoPhotonCoulomb} 
and~\ref{fig:Corr-bub-double}) enter at \text{N$^2$LO}\xspace.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.85\textwidth,clip]{TwoPhotonCoulomb}
\caption{Two-photon three-nucleon diagrams.  Solid lines represent nucleons 
whereas a deuteron (\ensuremath{{}^1S_0}\xspace dibaryon) is drawn as a double (thick) line.  Wavy 
lines are Coulomb photons.}
\label{fig:TwoPhotonCoulomb}
\end{figure}
\begin{figure}[htbp]
\centering
\includegraphics[width=0.2\textwidth]{Corr-bub-double}
\caption{Two-photon bubble.  Symbols are as in Fig.~\ref{fig:TwoPhotonCoulomb}.}
\label{fig:Corr-bub-double}
\end{figure}

The basic ingredient for all these topologies is the two-photon loop diagram 
shown in Fig.~\ref{fig:TwoPhotonExchange}.  This could be calculated directly, 
but it is more convenient to simply extract it as the $\mathcal{O}(\alpha^2)$ piece of 
the full off-shell Coulomb T-matrix.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.2\textwidth]{TwoPhotonExchange}
\caption{General two-photon exchange diagram.  Symbols are as in 
Fig.~\ref{fig:TwoPhotonCoulomb}.}
\label{fig:TwoPhotonExchange}
\end{figure}

For general kinematics with incoming (outgoing) momentum $\mathbf{p}$ ($\mathbf{q}$) and 
center-of-mass energy $E = k^2M_N$, this T-matrix can be written in the Hostler 
form~\cite{Chen:1972ab,Hostler:1964aa,Hostler:1964ab}
\begin{equation}
 T_C(k;\mathbf{p},\mathbf{q}) = \frac{\!4\pi\alpha}{(\mathbf{p}-\mathbf{q})^2}
 \left\{1-2\mathrm{i}\eta
 \int_1^\infty\left(\frac{s+1}{s-1}\right)^{\!\!-\mathrm{i}\eta}
 \frac{\mathrm{d} s}{s^2-1-\epsilon} \right\}
 \mathtext{,}
 \epsilon = \frac{(p^2-k^2)(q^2-k^2)}{k^2(\mathbf{p}-\mathbf{q})^2}
 \,,
\label{eq:TC-int}
\end{equation}
where
\begin{equation}
 \eta(k) = \frac{\alphaM_N}{2k}
\label{eq:eta-k}
\end{equation}
This form is most suitable for extracting the two-photon piece because it is 
straightforward to expand in $\alpha$.  Noting that $\eta=\mathcal{O}(\alpha)$ gives
\begin{equation}
 T_C^{(2)}(k;\mathbf{p},\mathbf{q}) = ({-}2\mathrm{i}\eta)\frac{4\pi\alpha}{(\mathbf{p}-\mathbf{q})^2}
 \int_1^\infty \frac{\mathrm{d} s}{s^2-1-\epsilon}
 = {-}\mathrm{i} \frac{4\pi\alpha^2M_N}{k(\mathbf{p}-\mathbf{q})^2}
 \frac{\mathrm{atanh}\left(\sqrt{1+\epsilon}\right)}{\sqrt{1+\epsilon}} \,.
\label{eq:TC2}
\end{equation}
A numerical comparison with the bare expression obtained from 
Fig.~\ref{fig:TwoPhotonExchange}---which involves a brute-force 
three-dimensional momentum integration---gives excellent agreement with 
Eq.~\eqref{eq:TC2}.

The diagrams shown in Fig.~\ref{fig:TwoPhotonCoulomb} can now be expressed in 
terms of $T_C^{(2)}$, based on previous results including the full Coulomb 
T-matrix~\cite{Konig:2014ufa} (\textit{cf.}\xspace~also 
Refs.~\cite{Kok:1979aa,Kok:1981aa,Ando:2010wq}).  Explicit expressions are 
given in Appendix~\ref{sec:CoulombKernels}.

Only the diagram with two photons exchange \emph{inside} a bubble, shown in 
Fig.~\ref{fig:Corr-bub-double}, requires some further work.  In principle, it 
can be obtained by integrating $T_C^{(2)}(k;\mathbf{p},\mathbf{q})$ over $\mathbf{p}$ and $\mathbf{q}$, 
but that would be unnecessarily tedious.  More conveniently, it can be 
extracted from the full Coulomb Green's 
function~\cite{Hostler:1963zz,Meixner:1933aa},
\begin{equation}
 G_C(E;\mathbf{r}'=\mathbf{0},\mathbf{r})
 = {-}2\mu\frac{\Gamma(1-\mathrm{i}\eta)\,W_{\mathrm{i}\eta;1/2}(-2\mathrm{i} k r)}{4\pi r}
 \mathtext{,} \eta = \frac{\alpha\mu}{k}
 \mathtext{,} \mu = \frac{M_N}{2} \,.
\label{eq:GC-W}
\end{equation}
in the limit $\mathbf{r}\to0$.  When this is expanded in $\alpha$, the first to 
terms are divergent, but the $\mathcal{O}(\alpha^2)$ and all higher orders are 
finite~\cite{Kong:1999sf,Kong:2000px}.  Expanding Eq.~\eqref{eq:GC-W} first in 
$r$ and subsequently in $\alpha$ gives\footnote{The logarithmic piece here 
actually comes together with a $\log r$ that has the same prefactor, and should 
thus be disregarded.}
\begin{equation}
 \frac{1}{2\mu} G_C(E;\mathbf{r}'=\mathbf{0},\mathbf{r}) = \text{div.\ terms}
 - \frac{\mathrm{i} k}{4\pi}
 + \frac{\alpha\mu\big[\log({-}2\mathrm{i} k) + C_E - 1\big]}{2\pi}
 - \frac{\mathrm{i} \mu^2 \pi \alpha^2}{12 k}
 + \mathcal{O}(r,\alpha^3) \,,
\label{eq:GC-W-expansion}
\end{equation}
such that one can read off
\begin{equation}
 \delta J_0^{(2)}(k) = {-}\frac{M_N}{4\pi}
 \left[\frac{\mathrm{i}\pi^2 \alpha^2M_N^2}{12k}\right] \,.
\label{eq:delta-J0-2}
\end{equation}
The label $\delta J_0^{(2)}(k)$ was chosen because alternatively this piece
can be extracted from the Coulomb bubble summing up two and more 
photon exchanges, which Ref.~\cite{Konig:2015aka} showed to be
\begin{equation}
 \delta J_0(k) = {-}\frac{\alphaM_N^2}{4\pi}
 \left[\psi(\mathrm{i}\eta) + \frac{1}{2\mathrm{i}\eta} + C_\Delta\right]
 + \frac{M_N}{4\pi} \mathrm{i} k \,.
\label{eq:delta-J0}
\end{equation}
Expanding the digamma function as
\begin{equation}
 \psi(x) = {-}\frac{1}{x} - C_E + \frac{\pi^2x}{6} + \mathcal{O}(x^2)
\end{equation}
gives again Eq.~\eqref{eq:delta-J0} and provides a consistency check.  In fact, 
all higher-order Coulomb contributions can be obtained in the same way as 
described here.

\subsection{Propagator corrections}
\label{sec:TwoBodyParamaters}

Next-to-leading order corrections to the dibaryon propagators are have been 
discussed in great detail in Ref.~\cite{Konig:2015aka}.  They are exactly the 
same here for the expansion that only takes the \ensuremath{{}^1S_0}\xspace channel in the 
unitarity leading order.  For the ``full unitarity'' expansion introduced in 
Ref.~\cite{Konig:2016utl}, the \ensuremath{{}^3S_1}\xspace channel is taken in the same limit at 
\text{LO}\xspace, such that generically
\begin{subalign}[eq:Delta-dt-0-1]
 \mathrm{i}\Delta_{d,t}^{(0)}(p_0,\mathbf{p})
 &= \frac{-\mathrm{i}}{\sigma_{d,t}^{(0)} + y_{d,t}^2 I_0(p_0,\mathbf{p})} \\
 \mathrm{i}\Delta_{t}^{(1)}(p_0,\mathbf{p})
 &= \mathrm{i}\Delta_{d,t}^{(0)}(p_0,\mathbf{p})\times
 \bigg[{-}\mathrm{i}\sigma_{d,t}^{(1)}
 - \mathrm{i} c_{d,t}^{(1)} \left(p_0-\frac{\mathbf{p}^2}{4M_N}\right) \bigg]
 \times \mathrm{i}\Delta_{d,t}^{(0)}(p_0,\mathbf{p}) \,.
\end{subalign}
with
\begin{equation}
 \sigma_{d,t}^{(0)} = \frac{2\Lambda}{\pi}
 \mathtext{,}
 \sigma_{d,t}^{(1)} = {-}\frac{1}{a_{d,t}}
 \mathtext{,}
 c_{d,t}^{(1)} = \frac{M_N r_{d,t}}{2} \,.
\label{eq:sigma-dt-unitarity}
\end{equation}
In particular, while in the standard and \ensuremath{{}^1S_0}\xspace-unitarity \text{LO}\xspace the \ensuremath{{}^3S_1}\xspace 
propagator is matched to the effective range expansion around the deuteron 
pole, giving\footnote{The deuteron binding momentum is taken here as 
$\gamma_d=45.701~\ensuremath{\mathrm{MeV}}$~\cite{vanderLeun:1982aa} and the effective range is 
$\rho_d=1.765~\ensuremath{\mathrm{fm}}$~\cite{deSwart:1995ui}.}
\begin{equation}
 \sigma_d^{(1)} = \frac{\rho_d\gamma_d^2}{2} \mathtext{,}
 c_d^{(1)} = \frac{M_N\rho_d}{2} \,,
\label{eq:sigmad1}
\end{equation}
the ordinary expansion around zero energy is used in the full-unitarity 
case to treat both channels fully equivalently.  Either way, at \text{N$^2$LO}\xspace there are 
quadratic insertions of $\sigma_{d,t}^{(1)}$ and $c_{d,t}^{(1)}$,
\begin{equation}
 \mathrm{i}\Delta_{t}^{(2)}(p_0,\mathbf{p})
 = \mathrm{i}\Delta_{t}^{(0)}(p_0,\mathbf{p})\times
 \Bigg\{\bigg[{-}\mathrm{i}\sigma_{t}^{(1)}
 - \mathrm{i} c_{d,t}^{(1)} \left(p_0-\frac{\mathbf{p}^2}{4M_N}\right)\Bigg\}^{\!2} \,,
\label{eq:Delta-t-2}
\end{equation}
and of course this matches directly onto the expansion of the renormalized 
amplitude when the expressions from Eq.~\eqref{eq:sigma-dt-unitarity} are 
inserted.  Note that the generic dibaryon parameter $\sigma$ (in a given 
channel) is related the corresponding ``ordinary'' $C_0$ used in 
Sec.~\ref{sec:TwoBody-NNLO} via $\sigma = {-}{4\pi}/{(M_N C_0)}$.  This 
means that their formal expansions are different, and in general, for a fixed 
scattering length included at \text{NLO}\xspace, one has $\sigma_{t(,pp)}^{(n)} = 0$ for $n 
\geq 2$.\footnote{The same is true for the \ensuremath{{}^3S_1}\xspace channel expanded around 
zero momentum, but with the ERE around the deuteron pole, there are 
contributions to $\sigma_d$ from all orders.}  In the $nd$ sector, however, 
there is a correction $\sigma_{t,nn}^{(2)} = {-}1/\ensuremath{a_{t,nn}} + 1/a_t$ which shifts 
the inverse scattering length away from its isospin-symmetric value.

\medskip
In the $pp$ channel, there are Coulomb contributions as well.  Most notably,
the bubble with a single photon exchanged inside---denoted as $\delta I_0(k)$---
is logarithmically divergent, and this divergence is absorbed by 
setting~\cite{Konig:2015aka}
\begin{equation}
 \sigma_{t,pp}^{(1)} = {-}\frac{1}{a_C}
 - \alphaM_N\left(\log\dfrac{2\Lambda}{\alphaM_N}
 - C_\zeta - C_\Delta\right)
\label{eq:renorm-sigma-t-pp-1}
\end{equation}
with known constants $C_\zeta$ and $C_\Delta$.\footnote{As discussed in detail 
in Ref.~\cite{Konig:2015aka}, the term $C_\zeta$ cancels exactly against the 
contribution from $\delta I_0(k)$, leaving only $C_\Delta\approx0.579$ in the 
final expressions.}  \text{NLO}\xspace also includes isospin-symmetric range corrections 
given by $c_t^{(1)}$.  At \text{N$^2$LO}\xspace, there are quadratic insertions of 
$\sigma_{t,pp}^{(1)}$, $c_t^{(1)}$, and $\delta I_0(k)$, but in addition also 
the isospin-breaking range correction $c_{t,pp}^{(2)} = M_N(r_C-r_t)/2$ and the 
genuine two-photon contribution $\delta J_0^{(2)}(k)$ discussed in the previous 
section.  Overall, this amounts to
\begin{spliteq}
 \mathrm{i}\Delta_{t,pp}^{(2)}(p_0,\mathbf{p})
 &= \mathrm{i}\Delta_{t,pp}^{(0)}(p_0,\mathbf{p})\times
 \Bigg\{\bigg[{-}\mathrm{i}\sigma_{t,pp}^{(1)} 
 - \mathrm{i} c_t^{(1)} \left(p_0-\frac{\mathbf{p}^2}{4M_N}\right) \\
 &\hspace{9em}\null
 - \mathrm{i} y_t^2\delta I_0\!\left(\mathrm{i}\sqrt{\mathbf{p}^2/4-M_N p_0-\mathrm{i}\varepsilon}\right)\bigg]
 \times \mathrm{i}\Delta_{t,pp}^{(0)}(p_0,\mathbf{p})\Bigg\}^{\!2} \\
 &\null + \mathrm{i}\Delta_{t,pp}^{(0)}(p_0,\mathbf{p})\times\bigg[
  {-}\mathrm{i} c_{t,pp}^{(2)} \left(p_0-\frac{\mathbf{p}^2}{4M_N}\right) \\
 &\hspace{9em}\null
  - \mathrm{i}y_t^2\delta J_0^{(2)}\!\left(\mathrm{i}\sqrt{\mathbf{p}^2/4-M_N p_0-\mathrm{i}\varepsilon}\right)
 \bigg]\times \mathrm{i}\Delta_{t,pp}^{(0)}(p_0,\mathbf{p}) \,,
\label{eq:Delta-pp-N2LO-1}
\end{spliteq}
and inserting the various definitions this simplifies to
\begin{spliteq}
 \mathrm{i}\Delta_{t,pp}^{(2)}(p_0,\mathbf{p}) &= {-}\mathrm{i}\dfrac{
 \Bigg\{\dfrac{1}{a_C} - \alphaM_N\left[C_\Delta
 + \log\!\left(\dfrac{\alphaM_N}{2\tilde{k}(p_0,\mathbf{p})}\right)
 \right] + \dfrac{r_t}{2}\,\tilde{k}(p_0,\mathbf{p})^2 \Bigg\}^{\!2}
 - \dfrac{\pi^2\alpha^2M_N^2}{12}}
 {\tilde{k}(p_0,\mathbf{p})^3} \\
 &\null - \mathrm{i} \frac{r_C - r_t}{2} \,,
\label{eq:Delta-pp-N2LO-2}
\end{spliteq}
where $\tilde{k}(p_0,\mathbf{p}) \equiv \sqrt{{\mathbf{p}^2}/{4}-M_N p_0-\mathrm{i}\varepsilon}$.

\section{Three-body observables up to \text{N$^2$LO}\xspace}
\label{sec:ThreeBody}

\subsection{Integral equations}
\label{sec:ThreeBodyEquations}

Adapting the notation of Refs.~\cite{Konig:2014ufa,Konig:2015aka}, the LO 
amplitude is a three-vector in channel space (with the last two components 
corresponding to $np$ and $pp/nn$ singlet-dibaryon legs):
\begin{equation}
 \vec{\mathcal{T}}_\sss \equiv \left(\TS^\mathrm{d,a},\TS^\mathrm{d,b1},\TS^\mathrm{d,b2}\right)^T \,.
\end{equation}
It is determined by the integral equation
\begin{equation}
 \vec{\mathcal{T}}_\sss^{(0)}
 = {-}\hat{K}^{(0)}\vec{e}_1
 + (\hat{K}^{(0)}\hat{D}^{(0)}) \otimes \vec{\mathcal{T}}_\sss^{(0)} \,,
\label{eq:Nd-IntEq-0}
\end{equation}
where $\vec{e}_1=(1,0,0)^T$ and $\otimes$ represents an integral over the 
intermediate momentum,
\begin{equation}
 A \otimes B \equiv \frac1{2\pi^2}
 \int_0^\Lambda\mathrm{d} q\,q^2\,A(\ldots,q)B(q,\ldots) \,,
\label{eq:SH-int}
\end{equation}
involving a momentum cutoff $\Lambda$.  $\hat{D}^{(0)})$ is a diagonal matrix of 
dibaryon propagators, \textit{viz.}\xspace
\begin{equation}
 \hat{D}^{(n)} = \mathrm{diag}\left(D_d^{(n)},D_t^{(n)},D_{t(,pp)}^{(n)}\right) \,,
\label{eq:D-n}
\end{equation}
where generically $D_x^{(n)}(E;q) \equiv 
\Delta_x^{(n)}\big(E-q^2/(2M_N);q\big)$, and the 
kernel matrix is
\begin{equation}
 \hat{K}^{(0)} = 2\pi \begin{pmatrix}
 {-}\left(K_\sss+\frac{2H_{0,0}(\Lambda)}{\Lambda^2}\right) &
 \left(3K_\sss+\frac{2H_{0,0}(\Lambda)}{\Lambda^2}\right) &
 \left(3K_\sss+\frac{2H_{0,0}(\Lambda)}{\Lambda^2}\right)\\[0.5em]
 \left(K_\sss+\frac{2H_{0,0}(\Lambda)}{3\Lambda^2}\right) &
 \left(K_\sss-\frac{2H_{0,0}(\Lambda)}{3\Lambda^2}\right) &
 {-}\left(K_\sss+\frac{2H_{0,0}(\Lambda)}{3\Lambda^2}\right)\\[0.5em]
 \left(2K_\sss+\frac{4H_{0,0}(\Lambda)}{3\Lambda^2}\right) &
 {-}\left(2K_\sss+\frac{4H_{0,0}(\Lambda)}{3\Lambda^2}\right) &
 {-}\frac{4H_{0,0}(\Lambda)}{3\Lambda^2}
 \end{pmatrix} \,.
\label{eq:K-0}
\end{equation}
It describes the isospin-symmetric system without Coulomb contributions and 
thus only contains the one-nucleon exchange part $K_\sss$ and the \text{LO}\xspace 
three-nucleon force $H_{0,0}$.  For more detailed definitions, see 
Ref.~\cite{Konig:2014ufa}.

As indicated by the subscript ``$\mathrm{s}$,'' the \text{LO}\xspace equation only contains 
contributions from the strong interaction.  The additional inclusion of Coulomb 
effects perturbatively builds up the ``full'' amplitude.  Up to \text{N$^2$LO}\xspace, the
corresponding integral equations are
\begin{equation}
 \vec{\mathcal{T}}_\fff^{(1)}
 = {-}\hat{K}^{(1)}\vec{e}_1
 + \left[\hat{K}^{(1)}\hat{D}^{(0)} + \hat{K}^{(0)}\hat{D}^{(1)} \right]
  \otimes \vec{\mathcal{T}}_\fff^{(0)}
 + (\hat{K}^{(0)}\hat{D}^{(0)}) \otimes \vec{\mathcal{T}}_\fff^{(1)}
\label{eq:Nd-IntEq-1}
\end{equation}

and
\begin{multline}
 \vec{\mathcal{T}}_\fff^{(2)}
 = {-}\hat{K}^{(2)}\vec{e}_1
 + \left[\hat{K}^{(2)}\hat{D}^{(0)}
 + \hat{K}^{(1)}\hat{D}^{(1)} + \hat{K}^{(0)}\hat{D}^{(2)}
 \right] \otimes \vec{\mathcal{T}}_\fff^{(0)} \\
 + \left[\hat{K}^{(1)}\hat{D}^{(0)} + \hat{K}^{(0)}\hat{D}^{(1)} \right]
  \otimes \vec{\mathcal{T}}_\fff^{(1)}
 + (\hat{K}^{(0)}\hat{D}^{(0)}) \otimes \vec{\mathcal{T}}_\fff^{(2)} \,,
\label{eq:Nd-IntEq-2}
\end{multline}
with $\vec{\mathcal{T}}_\fff^{(0)} = \vec{\Tgen_\sss}^{(0)}$.  Essentially, the structure is the 
same as in the two-body example discussed in Sec.~\ref{sec:TwoBody-NNLO}, only 
that here there are two types of corrections:
\begin{enumerate}
\item Contributions to the dibaryon propagators (from finite scattering 
lengths and effective ranges as well as two-body Coulomb effects in the $pp$ 
subsystem), encoded in the general expression given in Eq.~\eqref{eq:D-n}.
\item  The associated three-body force corrections as well as three-body Coulomb 
diagrams contributing to the higher-order interaction-kernel matrices.
\end{enumerate}
Specifically, these matrices are
\begin{equation}
 \hat{K}^{(1)} = 2\pi \begin{pmatrix}
 {-}\left(K_\text{bub}+K_\text{box}+\frac{2H_{0,1}(\Lambda)}{\Lambda^2}\right) &
 \left(3K_\text{box}+\frac{2H_{0,1}(\Lambda)}{\Lambda^2}\right) &
 \left(3K_\text{tri}^{(\text{out})}+\frac{2H_{0,1}(\Lambda)}{\Lambda^2}\right)\\[0.5em]
 \left(K_\text{box}+\frac{2H_{0,1}(\Lambda)}{3\Lambda^2}\right) &
 {-}\left(K_\text{bub}-K_\text{box}+\frac{2H_{0,1}(\Lambda)}{3\Lambda^2}\right) &
 {-}\left(K_\text{tri}^{(\text{out})}+\frac{2H_{0,1}(\Lambda)}{3\Lambda^2}\right)\\[0.5em]
 \left(2K_\text{tri}^{(\text{in})}+\frac{4H_{0,1}(\Lambda)}{3\Lambda^2}\right) &
 {-}\left(2K_\text{tri}^{(\text{in})}+\frac{4H_{0,1}(\Lambda)}{3\Lambda^2}\right) &
 {-}\frac{4H_{0,1}(\Lambda)}{3\Lambda^2}
 \end{pmatrix}
\label{eq:K-1}
\end{equation}
and
\begin{multline}
 \hat{K}^{(2)} = 2\pi \begin{pmatrix}
 {-}\left(K_\text{bub}^{(2)}+K_\text{box}^{(2)}+2H^{(2)}\right) &
 \left(3K_\text{box}+2H^{(2)}\right) &
 \left(3K_\text{tri}^{(\text{2,out})}+2H^{(2)}\right)\\[0.5em]
 \left(K_\text{box}^{(2)}+\frac23H^{(2)}\right) &
 {-}\left(K_\text{bub}^{(2)}-K_\text{box}^{(2)}+\frac23H^{(2)}\right) &
 {-}\left(K_\text{tri}^{(\text{2,out})}+\frac23H^{(2)}\right)\\[0.5em]
 \left(2K_\text{tri}^{(\text{2,in})}+\frac43H^{(2)}\right) &
 {-}\left(2K_\text{tri}^{(\text{2,in})}+\frac43H^{(2)}\right) &
 {-}\frac43H^{(2)}
 \end{pmatrix} \\
 \null + 2\pi \, \mathrm{diag}\big(K_{\rho_d},K_{r_t},0\big) \,,
\label{eq:K-2}
\end{multline}
where the collection of \text{N$^2$LO}\xspace three-nucleon forces is
\begin{equation}
 H^{(2)}(\Lambda) = \frac{H_{0,2}(\Lambda)}{\Lambda^2}
 + \frac{H_{2,2}(\Lambda)}{\Lambda^4}(M_N E + \gamma_d^2)
 + \frac{H_{0,2}^{(\alpha)}(\Lambda)}{\Lambda^2} \,.
\end{equation}

The functions $K_\text{bub}^{(2)}$, $K_\text{box}^{(2)}$, and $K_\text{tri}^{(\text{2,out})}$ correspond, 
respectively, to the three diagrams shown in Fig.~\ref{fig:TwoPhotonCoulomb}.  
$K_\text{tri}^{(\text{2,in})}$ comes from the reversed version of 
Fig.~\ref{fig:TwoPhotonCoulomb}(c), where the photon exchanges are on the left 
side of the diagram.  Detailed expressions for all these functions are given in 
the Appendix~\ref{sec:CoulombKernels}.  The analogous one-photon functions in 
Eq.~\eqref{eq:K-1} are defined (for example) in Ref.~\cite{Konig:2015aka} and 
are thus not repeated here.  Finally, the contributions $K_{\rho_d}$ and $K_{r_t}$, 
corresponding to a direct coupling of photons to the dibaryons,
\begin{equation}
 K_{r}(E;k,p) = {-}\frac{\alphaM_N r}{2kp}
 \,Q_0\!\left(-\frac{k^2+p^2+\lambda^2}{2kp} \right)
 \mathtext{,}
 Q_0(x) = \frac{1}{2}\log\left(\frac{x+1}{x-1}\right)
 \mathtext{,}
 r = \rho_d, r_t \,,
\label{eq:K_r}
\end{equation}
enter at \text{N$^2$LO}\xspace in the power counting ($\lambda$ here is a small photon mass 
introduced as infrared regulator and discussed further in the next section).

Note that the new three-body force $H_{0,2}^{(\alpha)}(\Lambda)$ is needed 
due to terms mixing range corrections and $O(\alpha)$ Coulomb contributions.
While from the \text{N$^2$LO}\xspace calculation based on the T-matrix method it is not 
directly obvious, one can identify contributions to the energy shift of the form
\begin{equation}
 \vcenteredhbox{%
  \raisebox{0.2em}{\includegraphics[width=5em]{DeltaE-bubble-rd}}
 }
 \sim \; \left(\frac1q\right)^{\!2}_{\!\!\mathcal{B}}
 \times \Big(q^3\Big)^{3}_{\!\text{loops}}
 \times \left(\frac1{q^2}\right)^{\!2}_{\!\!N}
 \times \left(\frac1q\right)^{\!3}_{\!\!D}
 \times \left(\frac1{q^2}\right)_{\!\text{photon}}
 \times \Big(q^2\Big)_{\!\text{range corr.}}
 \;\sim\; q^0 \,,
\end{equation}
where all loop momenta have been generically denoted by $q$ and the notation of 
Ref.~\cite{Konig:2014ufa} has been adapted to count ultraviolet behavior 
diagram.  This, as well as similar diagrams involving other topologies and/or 
a direct dibaryon-photon coupling, are logarithmically divergent, such that the 
three-body forces determined from the $nd$ system alone are no longer sufficient 
to ensure renormalization.

\subsection{Extraction of binding energies}

Binding energies can be obtained by plugging the solutions of the integral 
equations into the expressions given in Sec.~\ref{sec:SecondOrderEnergies}.
At LO, this merely amounts to varying the energy until the homogeneous version 
of Eq.~\eqref{eq:Nd-IntEq-0} has a solution.  When the \text{LO}\xspace three-nucleon force 
is fitted to the physical triton binding energy, the energy is fixed at that 
value and $H_{0,0}(\Lambda)$ is varied (for a given value of $\Lambda$).

The perturbative corrections to the binding energy are defined in terms of 
limits where the energy approaches the \text{LO}\xspace pole.  Numerically, this is evaluated 
by varying the energy around the pole and interpolating with a low-order 
polynomial to extract the residues (a linear interpolation is found to be 
sufficient in most cases if a 1\% interval around the pole is considered).
The reliability of this procedure has been checked by comparing \text{NLO}\xspace results 
with those obtained as matrix elements between leading-order vertex functions 
(as described, for example, in Ref.~\cite{Konig:2015aka}), finding excellent 
agreement.

\subsection{Subtracted phase shifts}

For $pd$ scattering it is necessary to calculate Coulomb-subtracted phase 
shifts $\tilde{\delta}(k) \equiv \delta_\mathrm{full}(k) - \delta_\mathrm{c}(k)$, where the 
first term---corresponding to the combination of strong and Coulomb 
contributions---is extracted perturbatively from the full amplitude 
(specifically, from the the $pd \to pd$ component $\Tgen_\fff\equiv\TF^\mathrm{d,a}$).  Up to 
\text{N$^2$LO}\xspace, the expressions are
\begin{subalign}
 \delta_{\mathrm{full}}^{(0)}(k)
 &= \frac{1}{2\mathrm{i}} \log\!\left(1+\tfrac{\mathrm{i} k\mu}{\pi}
 T_\mathrm{full}^{(0)}(k)\right) \,, \\
 \delta_{\mathrm{full}}^{(1)}(k) &= \frac{k\mu}{2\pi}
 \dfrac{T_\mathrm{full}^{(1)}(k)}{1+\tfrac{\mathrm{i} k\mu}{\pi} T_\mathrm{full}^{(0)}(k)} \,, \\
 \delta_{\mathrm{full}}^{(2)}(k) &= \frac{k\mu}{2\pi}
 \dfrac{T_\mathrm{full}^{(2)}(k)}{1+\tfrac{\mathrm{i} k\mu}{\pi} T_\mathrm{full}^{(0)}(k)}
 - \mathrm{i} \delta_{\mathrm{full}}^{(1)}(k)^2 \,,
\end{subalign}
where the T-matrices, in turn, combine the perturbative expansion of the 
deuteron wavefunction renormalization,
\begin{equation}
 Z_0^{-1} = \mathrm{i}\frac{\partial}{\partial p_0}
 \left.\frac{1}{\Delta_d(p_0,\mathbf{p})}\right|_{p_0 =-\frac{\gamma_d^2}{M_N},\,\mathbf{p}=0}
 \implies Z_0 = Z_0^{(0)} + Z_0^{(1)} + \cdots \,,
\label{eq:Z0}
\end{equation}
with that of the bare amplitude, \textit{viz.}\xspace
\begin{equation}
 T_\mathrm{full}^{(0)}(k) = Z_0^{(0)}\Tgen_\fff^{(0)}(E_k;k,k)
 \mathtext{,}
 T_\mathrm{full}^{(1)}(k) = Z_0^{(0)}\Tgen_\fff^{(1)}(E_k;k,k)+Z_0^{(1)}\Tgen_\fff^{(0)}(E_k;k,k) 
\,,
\end{equation}
and following the same pattern for $\Tgen_\fff^{(2)}(k)$.  As before, the reduced mass 
is $\mu=2M_N/3$.

Most Coulomb contributions are formally divergent in the on-shell limit.  This 
infrared singularity corresponds to the long-range nature of the Coulomb force 
and is conveniently regulated by introducing a small photon mass $\lambda$ 
where necessary, as it was done for example in Eq.~\eqref{eq:K_r}.  While in 
the low-energy limit, where Coulomb effects are highly nonperturbative, this 
``screening'' is very important~\cite{Konig:2013cia}, in the perturbative 
regimes considered here (\ensuremath{{}^3\mathrm{He}}\xspace and $pd$ scattering at intermediate 
center-of-mass momenta), there is almost no noticeable $\lambda$ 
dependence in the subtracted phase shifts for photon masses below about 
$1~\ensuremath{\mathrm{MeV}}$ (\textit{cf.}\xspace~results shown in Fig.~\ref{fig:Phase-pd-U-2-lambda}). In practice,
a fixed small value can be used instead of numerically taking the limit
$\lambda\to0$.

According to Ref.~\cite{Konig:2013cia}, the pure Coulomb part should be taken 
in a simple two-body picture in order to compare to experimental phase-shift 
analyses and potential-model calculations.   Including the photon mass, this 
corresponds to a Yukawa interaction between proton and deuteron.  At \text{LO}\xspace, there 
are no Coulomb effects by construction, so $T_\mathrm{c}^{(0)}(k)=0$ and 
$\delta_{\mathrm{c}}^{(0)}(k) = 0$.  At \text{NLO}\xspace and \text{N$^2$LO}\xspace, exact expressions can be 
obtained from the analogous case of pion-exchange contributions in the \ensuremath{{}^1S_0}\xspace 
nucleon-nucleon channel.  Adapting results obtained 
in Ref.~\cite{Fleming:1999ee} for pion-exchange diagrams to the 
formally analogous $pd$-Coulomb case gives
\begin{equation}
 \delta_\mathrm{c}^{(1)}(k) = \frac{k\mu}{2\pi} T_\mathrm{c}^{(1)}(k)
 = {-}\frac{\alpha\mu}{2k}\log\left(1+\frac{4k^2}{\lambda^2}\right)
\end{equation}
and $\delta_\mathrm{c}^{(2)}(k) = k\mu\,T_\mathrm{c}^{(2)}(k)/(2\pi) - 
\mathrm{i}\delta_\mathrm{c}^{(1)}(k)^2$
with
\begin{equation}
 T_\mathrm{c}^{(2)}(k) = \frac{2\mu\alpha^2\pi}{k^3}\left[
 \frac{\mathrm{i}}{4}\log\left(1+\frac{4k^2}{\lambda^2}\right)^{\!2}
 + \mathrm{Im}\,\mathrm{Li}_2\!\left(\frac{2k^2 - \mathrm{i} k\lambda}{\lambda^2+4k^2}\right)
 + \mathrm{Im}\,\mathrm{Li}_2\!\left(\frac{{-}2k^2 + \mathrm{i} k\lambda}{\lambda^2}\right)
 \right] \,.
\end{equation}

\section{Results}
\label{sec:Results}

\subsection{The $\boldsymbol{n}\boldsymbol{d}$ system}

As a first step, the \ensuremath{{}^1S_0}\xspace-unitarity expansion can be tested at \text{N$^2$LO}\xspace in the 
$nd$ system.  The corresponding amplitudes are readily obtained by omitting all 
Coulomb contributions in the integral equations discussed in 
Sec.~\ref{sec:ThreeBodyEquations} and changing the $pp$ to an $nn$ propagator.  
The $nd$ doublet scattering length can then be calculated as
\begin{equation}
 \ensuremath{{}^2a_{\text{$n$--$d$}}} = {-}\frac{M_N}{3\pi}\lim\nolimits_{k\to0} Z_0\TS^\mathrm{d,a}(E_k;k,k) \,,
\label{eq:and-2}
\end{equation}
with appropriate perturbative expressions at each order.

While Ref.~\cite{Konig:2015aka} found that already at \text{NLO}\xspace in the $1/|a_t|$ 
expansion, results for both the Phillips line (\textit{i.e.}, the correlation between 
$\ensuremath{{}^2a_{\text{$n$--$d$}}}$ and the \Triton binding energy~\cite{Phillips:1968zze}) as well 
as $nd$ phase shifts fall almost on top of those in the standard pionless 
counting (including $1/a_t$ nonperturbatively at \text{LO}\xspace), the agreement is even 
better at \text{N$^2$LO}\xspace (relative differences compared to standard \text{LO}\xspace are an 
order of magnitude smaller at \text{N$^2$LO}\xspace than at \text{NLO}\xspace, of the order $10^{-5}$ for 
the Phillips line and $10^{-4}$ for the $nd$ phase shift), thus establishing 
more firmly that the expansion converges rapidly.

\subsection{${}^{\boldsymbol{3}}\mathbf{He}$ binding energy}
\label{sec:ThreeHeResults}

Studying the \ensuremath{{}^3\mathrm{He}}\xspace bound state is more interesting.  Setting all range 
corrections to zero, it is possible to compare the expansions suggested, 
respectively, in Refs.~\cite{Konig:2015aka} and~\cite{Konig:2016utl} at \text{N$^2$LO}\xspace.  
In the first case, only the \ensuremath{{}^1S_0}\xspace channel is expanded around the unitarity 
limit, whereas the latter does the same for the \ensuremath{{}^3S_1}\xspace channel.  Since 
$a_d\sim5.4~\ensuremath{\mathrm{fm}} < |a_t|\sim23.7~\ensuremath{\mathrm{fm}}$, corrections are expected to be larger in 
the full-unitarity case.

\begin{figure}[tb]
\centering
\includegraphics[width=0.75\textwidth,clip]{En-3He-UU-2}
\caption{%
Three-nucleon binding energies up to second-order in the scattering-length 
expansion (\textit{i.e.}, without effective-range corrections).  The (red) dashed curve 
and diamonds show results for the expansion of Ref.~\cite{Konig:2015aka}, which 
takes only the \ensuremath{{}^1S_0}\xspace channel in the unitarity limit at \text{LO}\xspace.  The (green) 
dotted line and circles use the ``full-unitarity'' scheme of 
Ref.~\cite{Konig:2016utl}, which takes the \ensuremath{{}^3S_1}\xspace channel in the 
\text{LO}\xspace-unitarity limit as well.  For each calculation, the isospin-symmetric 
three-nucleon forces have been fitted to keep the triton binding energy fixed 
at its physical value.  The horizontal dotted line indicates the experimental 
value for the \ensuremath{{}^3\mathrm{He}}\xspace binding energy.}
\label{fig:En-3He-UU-2}
\end{figure}

Results are shown in Fig.~\ref{fig:En-3He-UU-2}.  Note that convergence towards 
the experimental value is not to be expected in this expansion that (besides 
Coulomb effects) only perturbatively includes the finite scattering length(s).  
What is important, however, is that the results for full and \ensuremath{{}^1S_0}\xspace-only 
unitarity move closer together at \text{N$^2$LO}\xspace, which indicates that both 
expansions---and in particular the perturbative treatment of $1/a_d$---are 
indeed convergent.

While for \ensuremath{{}^1S_0}\xspace unitarity there is a little less binding at \text{N$^2$LO}\xspace compared 
to \text{NLO}\xspace, the effect is opposite for the full-unitarity \text{LO}\xspace case.  This can 
generally be understood as a cancellations between different contributions at 
this order, which play out slightly differently in the two cases.  Starting at 
\text{NLO}\xspace in the full-unitarity scheme, $1/a_d$ corrections add attraction (in the 
unitarity limit, there is less attraction in the \ensuremath{{}^3S_1}\xspace channel compared to 
finite $a_d>0$), requiring a repulsive contribution to the three-nucleon force 
to keep the triton in place.  The opposite effect happens in the \ensuremath{{}^1S_0}\xspace 
channel (where the physical scattering length is negative), but it is smaller 
compared to the triplet contribution since $1/|a_t| < 1/a_d$.  Coulomb effects 
and $1/a_C<0$ provide pure repulsion at \text{NLO}\xspace, and this is the only net effect 
at that order because all isospin-symmetric contributions cancel when the 
\Triton--\ensuremath{{}^3\mathrm{He}}\xspace binding-energy difference is calculated in first-order 
perturbation theory.  At \text{N$^2$LO}\xspace there are nonlinear (quadratic) mixtures of the 
various \text{NLO}\xspace insertions, as well as genuine \text{N$^2$LO}\xspace terms (two-photon diagrams 
and the \ensuremath{a_{t,nn}} correction) that enter linearly at this order, so that there is no 
clear \textit{a priori}\xspace intuition what the net effect should be.  In any case, the 
overall corrections are small (note the zoomed scale in 
Fig.~\ref{fig:En-3He-UU-2}), as one would expect at \text{N$^2$LO}\xspace in this expansion.

\medskip
Range corrections are generally an important effect, but they cancel at \text{NLO}\xspace 
for the same reason as discussed above (adjustment of the three-nucleon force). 
At \text{N$^2$LO}\xspace, their inclusion gives rise to divergences, stemming both directly 
from isospin breaking in the effective ranges, $r_C \neq r_t$ as well as from 
the interference of (isospin-symmetric) range corrections and Coulomb 
contributions that entered together at \text{NLO}\xspace.  An isospin breaking three-nucleon 
force is thus required at this order for renormalization (just like the \text{NLO}\xspace 
term identified in Ref.~\cite{Vanasse:2014kxa} for standard pionless counting 
with Coulomb effects included nonperturbatively), and it is no longer possible 
to determine the \ensuremath{{}^3\mathrm{He}}\xspace binding energy without input from a $pd$ doublet 
observable.  This input can be chosen to be the \ensuremath{{}^3\mathrm{He}}\xspace energy itself, 
allowing a prediction of $pd$ scattering.

\subsection{$\boldsymbol{p}\boldsymbol{d}$ scattering}
\label{sec:pdResults}

Scattering phase shifts for the $pd$ doublet channel are shown in 
Fig.~\ref{fig:Phase-pd-U}.  For this calculation, the three-nucleon 
force fitting strategy has been adjusted to match what was used in 
Ref.~\cite{Vanasse:2013sda}, \textit{i.e.}, instead of determining $H_{0,n}(\Lambda)$ by 
demanding that the triton remains at its physical energy, the $nd$ doublet 
scattering length---calculated according to Eq.~\eqref{eq:and-2}---is used as 
input instead.  The triton is then moved into place at \text{N$^2$LO}\xspace by adjusting 
$H_{2,2}(\Lambda)$.\footnote{At \text{NLO}\xspace, it is also possible to calculate \ensuremath{{}^3\mathrm{He}}\xspace 
with $H_{0,n}(\Lambda)$ fitted to the $nd$ scattering length.  While this 
would move the absolute position of the predicted binding energy, the 
difference to the triton---which is what is really calculated---remains 
unchanged.}  On top of this, $H_{0,2}^{(\alpha)}(\Lambda)$ is finally adjusted 
to reproduce the \ensuremath{{}^3\mathrm{He}}\xspace energy.  Following this procedure has the advantage 
that all terms can be determined independently in subsequent one-parameter fits. 
 Note that this calculation is carried out for the \ensuremath{{}^1S_0}\xspace-unitarity scheme 
only in order to have a fixed deuteron bound state.

\begin{figure}[tb]
\centering
\includegraphics[width=0.75\textwidth,clip]{Phase-pd-U}
\caption{%
Proton-deuteron S-wave phase shifts calculated up to \text{N$^2$LO}\xspace, with the \text{LO}\xspace 
calculation identical to the $nd$ result and Coulomb and other isospin breaking 
effects added perturbatively at subsequent orders.  The isospin-breaking 
three-nucleon force has been fitted at to reproduce the experimental value for 
the \ensuremath{{}^3\mathrm{He}}\xspace binding energy at \text{N$^2$LO}\xspace.  All curves were calculated at a cutoff 
$\Lambda=4.8~\ensuremath{\mathrm{GeV}}$ and a fixed photon mass $\lambda=0.25~\ensuremath{\mathrm{MeV}}$ has been used as 
infrared regulator.  As demonstrated in Fig.~\ref{fig:Phase-pd-U-2-lambda}, 
which shows the dependence of the \text{N$^2$LO}\xspace result, this is sufficiently close to the 
limit $\lambda\to0$ for practical purposes.  The filled circles were calculated 
in Ref.~\cite{Kievsky:1996ca} using the AV18+UIX potential model.  Experimental 
data (open diamonds) is shown from Ref.~\cite{Arvieux:1974fma}.}
\label{fig:Phase-pd-U}
\end{figure}
\begin{figure}[tb]
\centering
\includegraphics[width=0.75\textwidth,clip]{Phase-pd-U-2-lambda}
\caption{%
Proton-deuteron S-wave phase shifts calculated at \text{N$^2$LO}\xspace for different values of 
the regulating photon mass $\lambda$.  The logarithm is taken to enhance 
the visibility of small differences, and overall convergence as $\lambda$ is 
lowered towards zero is apparent.  All curves were calculated at a cutoff 
$\Lambda=4.8~\ensuremath{\mathrm{GeV}}$}
\label{fig:Phase-pd-U-2-lambda}
\end{figure}

The \text{LO}\xspace result (green dotted curve in Fig.~\ref{fig:Phase-pd-U}) coincides 
with that for $nd$ scattering.  Since all Coulomb effects are included 
perturbatively here, the $pd$ calculation is not expected to work at very low 
center-of-mass momenta $k$, but already at $k=20~\ensuremath{\mathrm{MeV}}$ the Coulomb parameter 
$\alphaM_N/k$ is about $1/3$ (in fact, for the $pd$ system one should rather 
count $\alpha2M_N/(3k)$ to account for the proper reduced mass, and this is 
$\simeq1/4$ at $k=20~\ensuremath{\mathrm{MeV}}$).  Consequently, the perturbative calculation should 
work well for $k\gsim20~\ensuremath{\mathrm{MeV}}$, and higher-order Coulomb effects are certainly 
not the primary source of uncertainty at \text{N$^2$LO}\xspace.  Following 
Ref.~\cite{Konig:2015aka}, the $Q/\Lambda_{\slashed\pi}$ and $\aleph_0/Q$ expansions are 
paired by taking $Q\sim\sqrt{\mathstrut\aleph_0\Lambda_{\slashed\pi}}$.  Conservatively 
choosing $\aleph_0 \sim 1/a_C$ (\textit{i.e.}, assuming it is dominated by the largest 
scale) and $\Lambda_{\slashed\pi} \sim M_\pi$ gives $Q\sim60~\ensuremath{\mathrm{MeV}}$, slightly above the 
standard pionless counting which generally assumes $Q\sim\gamma_d\sim50~\ensuremath{\mathrm{MeV}}$.  
With this, the uncertainties at \text{NLO}\xspace and \text{N$^2$LO}\xspace are estimated, respectively, as 
$Q((\aleph_0/Q)^2) \simeq 20\%$ and $Q((\aleph_0/Q)^3) \simeq 8\%$.  They are 
indicated as shaded bands in Fig.~\ref{fig:Phase-pd-U}.  Within these bands, 
potential-model results from Ref.~\cite{Kievsky:1996ca} (using the AV18+UIX 
interaction) are accommodated.  Results from an experimental phase-shift 
analysis~\cite{Arvieux:1974fma} are shown as well.  While overall the agreement 
looks reasonable, a more rigorous statement cannot be made in the absence of 
error bars for the data points.

\section{Summary and outlook}
\label{sec:Conclusion}

This work demonstrates that the pionless counting schemes developed in 
Refs.~\cite{Konig:2015aka} and~\cite{Konig:2016utl}, which take, respectively, 
the \ensuremath{{}^1S_0}\xspace or both $NN$ S-wave channels in the unitarity limit at leading 
order, work well up to second-order perturbation theory.  The inclusion of 
effective-range corrections requires an isospin-breaking three-nucleon force in 
the doublet channel at \text{N$^2$LO}\xspace, but fixing this to a single datum (like the 
\ensuremath{{}^3\mathrm{He}}\xspace binding energy) allows predictions to be made at other energies.  In 
particular, there exists a regime where $pd$ scattering is well described using 
only perturbative Coulomb effects (with two-photon contributions included at 
\text{N$^2$LO}\xspace).  Within the assigned \textit{a priori}\xspace uncertainty based on the EFT expansion, 
results agree with potential-model calculations of $pd$ scattering.  This 
establishes that \text{N$^2$LO}\xspace calculations of the nuclear three-body system including 
Coulomb effects are feasible and well-behaved, which is an important proof of 
principle on the way to similar calculations to be carried out in chiral EFT.

It is clear that some nonperturbative treatment of Coulomb effects is 
needed as $k\to0$ in the three-nucleon system, \textit{i.e.}, for a calculation of the 
$pd$ scattering length.  At the very least, the diagram with a single photon 
attached to the bubble should be iterated at \text{LO}\xspace to reflect its infrared 
enhancement at small energies.  For the quartet channel (total spin 
$\nicefrac32$), exactly this was done in Ref.~\cite{Konig:2013cia}, with the box 
diagram and the dibaryon-photon range coupling added perturbatively (note 
that treating Coulomb diagrams differently based on which momentum 
scale dominates their behavior is generally what has been suggested in 
Ref.~\cite{Rupak:2001ci}).  Extending this procedure to the doublet channel, 
using the \ensuremath{{}^1S_0}\xspace-unitarity expansion to handle the Coulomb divergence in the 
$pp$ channel, is left for future work, as is an analysis of Coulomb effects in 
$pd$ scattering above the deuteron breakup threshold ($k\simeq53~\ensuremath{\mathrm{MeV}}$), where 
the energy in the two-body subsystem goes through zero.  Moreover, it will be 
interesting for future work to combine the unitarity-\text{LO}\xspace expansions employed 
here with the method developed in Ref.~\cite{Vanasse:2015fph} in order to 
calculate the \ensuremath{{}^3\mathrm{He}}\xspace charge radius.

Typical for pionless EFT, the calculation carried out here allows for its 
various low-energy constants to be determined independently.  This is done
analytically in the two-nucleon sector and using subsequent one-parameter fits 
for the $3N$ forces, which is conveniently implemented numerically.  However, by 
not taking into account the global uncertainty at a given order, one effectively 
assumes some higher-order terms not to contribute.  In Ref.~\cite{Konig:2015aka} 
it was already noted that by fixing the $pp$ scattering length exactly to its 
experimental value at \text{NLO}\xspace one overestimates the role of Coulomb effects over 
short-range isospin-breaking contributions.  Something similar happens by 
demanding that the \ensuremath{{}^3\mathrm{He}}\xspace binding energy takes its exact physical value at 
\text{N$^2$LO}\xspace via the adjustment of the isospin-breaking three-nucleon force.  While of 
course such effects are reflected in the assigned uncertainty estimates, it 
would be interesting to carry out a more sophisticated fitting procedure that 
propagates the uncertainties and correlations through different orders and 
derives probability distributions for the low-energy constants.  Such methods 
are currently being developed, with the primary goal of applying them to chiral 
EFT (see for example Ref.~\cite{Wesolowski:2015fqa} and earlier references 
therein).  Pionless EFT provides a convenient test case for such a framework 
because results can be compared directly to explicitly determined values for the 
LECs.  Work along these lines is in progress.

\begin{acknowledgments}
I would like to thank R.~Furnstahl and U.~van Kolck for useful discussions and 
valuable comments on the manuscript.  Moreover, I am grateful to the 
participants of the GSI-funded EMMI RRTF workshop \emph{ER15-02: Systematic 
Treatment of the Coulomb Interaction in Few-Body Systems} for a very stimulating 
meeting.  This work is supported in part by the NSF under award PHY--1306250, 
by the DOE-funded NUCLEI SciDAC Collaboration under award DE-SC0008533, as well 
as by the ERC Grant No.\ 307986 STRONGINT.
\end{acknowledgments}

