\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}