\section{Introduction}
\label{sect:introduction}
A multi-messenger detector for cosmic rays, capable of detecting charged particles and gamma-ray photons, represents an appealing
perspective for the next generation of space-based, high-energy direct detection experiments from both the scientific and the logistic
points of view. Such a facility would provide invaluable data for various research fields, at a significantly reduced cost compared to a
classic two dedicated detector solution (one for gamma rays and one for charged particles). Its design however poses many challenges.
In order to make precise measurements with a good statistical significance also at high energies, scientific requirements range from good
energy and tracking resolution to high geometric factor and effective area. Trying to meet all these requirements for charged particles and
gamma rays with a single detector raises some inevitable conflicts.

Gamma-ray tracking is typically performed using a pair-conversion telescope. Each plane of the telescope is made of a layer of passive
material (e.g. tungsten) in which the primary photon is converted into an e$^+$-e$^-$ pair, coupled with a sensitive layer (e.g. silicon
microstrip detector) which measures the impact position of the electron and the positron. Several of these planes are stacked to increase
the photon conversion probability and to provide multiple position measurements along the e$^+$-e$^-$ trajectory; from these measurements
the track of the incident photon is reconstructed. The planes are suitably spaced apart from each other, so as to achieve the lever arm
that, in combination with the spatial resolution of the position measurements, provides the required tracking resolution. This concept has
been successfully implemented in recent space-borne gamma-ray detectors, like Fermi \cite{Atwood2007} and AGILE \cite{Bulgarelli2010}, and
can be considered as the reference technique for gamma-ray tracking at energies above some hundreds MeV. Being sensitive to ionization, such
a device is in principle also able to track charged primary particles and is actually used in detectors for charged particles with gamma-ray
detection capabilities like DAMPE \cite{Chang2017}

On the other hand, charged particle physics requires an accurate measurement of the charge of the nuclei. The
charge $Z$ of a relativistic particle or atomic nucleus impinging on a pair-conversion telescope can be measured by exploiting the $Z^2$
scaling of the energy deposited by ionization on the sensitive layers. However, the presence of a non-negligible amount of passive material
increases the fragmentation probability of the nucleus. After fragmentation, the subsequent ionization measurements are no longer useful for
charge determination, and typically must be rejected in offline charge reconstruction by an algorithm which identifies the interaction
point. The reduced set of independent ionization measurements and the limited accuracy of the algorithm likely lead to a worsening of
the charge identification performance. The simplest solution to this problem is to add a dedicated sub-detector for charge measurement
externally to the pair-conversion telescope, for example a set of silicon detectors with a minimal amount of passive material. In this way
the charge can be precisely measured before the nucleus traverses a considerable amount of material, at the price of an increase in power
consumption and cost.

Such a configuration consisting of a charge detector on top of a pair-conversion telescope can provide robust performance for both charge
identification and tracking of gamma rays and charged particles, but it will have an important impact on the overall geometric factor of
the instrument. The ability to resolve point sources with the desired accuracy requires a good tracking resolution, that translates into a
high lever arm value for the tracking telescope (since at fixed resolution of position measurement the angular resolution scales as
$1/L$, where $L$ is the distance between the outmost tracking layers); stacking the charge detector on top of the telescope further
increases the overall height of the instrument. This results in a reduction of the field of view and so in a degradation of the geometric
factor, since for a particle telescope with a distance $d$ between the entrance and exit planes defining the instrument acceptance the
geometric factor scales approximately as $1/d^2$ \cite{Sullivan1971}. This fact will impair the possibility to detect a sufficient amount
of high-energy particles, effectively limiting the energy range covered by the instrument.

The weight of the pair-conversion telescope can also have a detrimental impact on detector design. The mass budget of a space experiment is
usually severely constrained, limiting the maximum weight and dimension of the calorimeter and thus the overall geometric factor.
Allocating a significant fraction of this budget to the passive material needed for the conversion of the gamma rays further reduces the
instrument performance.

In this paper, an alternative detector design named Tracker-In-Calorimeter (TIC) aimed at mitigating the issues described above is
presented, together with expected performance. Section \ref{sect:measurementprinciple} describes the measurement principle; section
\ref{sect:detectordesign} gives an example of a realistic detector design implementing the measurement principle; section
\ref{sect:mcsimulations} shows the expected performance for a full-scale detector obtained with Monte Carlo simulations, while section
\ref{sect:prototype} illustrates the design of a detector prototype. The results obtained using data from the prototype acquired at test
beams are presented in section \ref{sect:tbresults}.

\section{Measurement principle}
\label{sect:measurementprinciple}

The basic idea is to use a calorimeter as a tracking device for gamma rays. It relies on the fact that most cosmic-ray detectors feature an
electromagnetic calorimeter below the main tracking device, with the main purpose of measuring the particle energy and discriminating
electrons and protons by means of topological shower analysis. Indeed, the effort towards the high-energy frontier in space (sub-TeV for
gamma rays, multi-TeV for electrons and PeV for hadrons) is being pursued by running \cite{Asaoka2017,Chang2017} and future
\cite{Zhang2017} experiments mainly with deep, homogenous electromagnetic calorimeters with scintillating crystals. When a photon hits
the calorimeter it generates an electromagnetic shower, whose axis preserves the information about its direction. By sampling the transverse
profile of the shower at different depths in the calorimeter it is possible to reconstruct the photon direction, e.g. by computing the
center-of-gravity of the shower image on each layer and then fitting with a line. The accuracy and precision of this method can be improved
by more sophisticated reconstruction algorithms and also by a more refined lateral sampling.

This tracking principle has already been implemented by other experiments, obtaining resolutions better than about 100 $\mu$m on the
reconstructed impact point for electrons above 100 GeV using a Si-W sampling calorimeter \cite{Adriani2010}. Since the track reconstruction
procedure relies completely on the development of an electromagnetic shower, similar performance can be expected for photons. A successful
application of this method to the calorimetric section of a cosmic-ray detector would allow to offload the gamma-ray tracking duties to the
calorimeter, optimising the tracker design for the tracking and possible charge measurement of primary nuclei and electrons. In practice, it
would be possible to remove the passive material from the main tracker and also to reduce the overall height of the instrument (since
tracking requirements for gamma rays are usually more demanding than those for charged particles due to the necessity to identify point
sources), thus providing a solution to most of the issues described in section \ref{sect:introduction}.

\section{Detector design}
\label{sect:detectordesign}
This proposal for a realistic configuration of a tracking calorimeter for a future multi-messenger cosmic-ray experiment starts by
considering the innovative proposal for a homogeneous and isotropic calorimeter made by the CaloCube collaboration
\cite{Adriani2017,Adriani2019}. The main feature of this calorimeter is to be able to accept particles not only from the face pointing
towards the zenith but also from lateral faces, thus boosting the geometric factor by a theoretical 5x factor in the same mass and power
consumption envelope. This design has been chosen with some modifications for the future HERD experiment \cite{Zhang2017}, and thus
constitutes an appealing starting point. The CaloCube sensitive elements are cubic scintillating crystal of size 1 Moli\`ere
radius read out by two photodiodes of different sizes in order to achieve a high dynamic range. The crystals are arranged in a
cubic 3D mesh. Such a segmented calorimeter already provides sizeable information about the development of the electromagnetic shower,
that can be further enhanced by instrumenting one side (e.g. the zenithal one) with silicon microstrip detectors interleaved with the first
layers of scintillating crystals, as schematically shown in figure \ref{fig:design}. Going from top to bottom, the upper layers are made of
thin crystals rather than cubic crystals, and a gap is left between the first two silicon detectors. This arrangement reduces the crystal
material traversed by the e$^+$-e$^-$ generated by photon conversion before hitting the precise silicon sensors, mitigating the impact of
multiple scattering on the tracking performance. In addition to this, the small gap allows for having two consecutive position measurements
by silicon detectors without material in between, further reducing the multiple scattering effects. The two subsequent layers of crystals
are made of thin elements as well, providing a limited amount of material for shower development between other precise position
measurements. Then the basic CaloCube design made of cubic scintillators is restored. The silicon detectors sample the lateral development
of the early stage of electromagnetic showers originating in the crystals at different depths in the calorimeter with a sampling interval of
the transverse positions much smaller than that of crystals (typically hundreds of $\mu$m against some cm). This fine-grained information is
crucial in improving the track reconstruction performance. A similar approach has been adopted in \cite{Adriani2010}, where silicon
detectors are interleaved with tungsten plates to build a sampling calorimeter with tracking capabilities. The novelty in this proposal is
the use of an active material as gamma-ray converter to maximize the energy resolution.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=13cm]{improveddesign.eps}
\end{center}
\vspace{-0.5cm}
\caption{Side view of a portion of the tracking calorimeter. The first 3 scintillator layers on top are made of thin crystals (blue)
and
interleaved with 4 silicon microstrip layers (orange). Cubic crystals (cyan) are used below this top structure.}
\label{fig:design}
\end{figure}

\section{Monte Carlo simulations}
\label{sect:mcsimulations}
To quantitatively assess the expected performance of the tracking calorimeter proposed in the previous section, a FLUKA
\cite{Battistoni2015,Bohlen2014} Monte Carlo simulation has been set up. The simulated detector is modeled adopting some of the design
choices of the future HERD calorimeter \cite{Zhang2017} (which in turn follows in many aspects the CaloCube design described in section
\ref{sect:detectordesign}): the thin and cubic crystals are made of Lutetium Yttrium Orthosilicate (Lu$_{2(1-x)}$Y$_{2x}$SiO$_5$), commonly
known as LYSO,
and the cubes have side 3 cm while the thin scintillators are 1.5 cm (i.e. half-cube) thick. The simulated geometry is implemented after
that depicted in figure \ref{fig:design}. The whole calorimeter is a 21x21x21 mesh of cubic crystals, with the silicon layers and thin
crystals stacked atop the top face. The gap $L$ between the first two silicon detectors is 7 cm, while the inter-crystal gaps are 5 mm. The
total depth along the up-down direction is $\sim$ 3.2 nuclear interaction lengths ($\lambda_I$) and $\sim$ 58 radiation lengths ($X_0$). The
thickness and readout pitch for silicon detectors are set to 300 $\mu$m and 50 $\mu$m, respectively. The gaps between the cubes are filled
with carbon fiber, to mimic the presence of a support structure; no other supports are present.

The simulated input spectra consist of monochromatic gamma rays at different energies (1, 10 and 100 GeV), with an isotropic angular
distribution and a zenith angle for the initial direction limited in the $\pm$ 30$^\circ$ region around the zenith direction. The starting
point is sampled from a flat surface just above the first layer of thin crystals, so that the whole top face is uniformly illuminated by
incoming photons.

\subsection{Track reconstruction algorithm}
\label{sect:algorithm}
The estimator of the incident particle track is an iterative algorithm based on all the information obtained by the calorimeter signals
(both crystals and Si detectors). The key point is that the energy of the primary photon is well known thanks to the calorimetric
measurement, therefore all the parameters of the algorithm, which to a certain extent can depend on the energy, can be considered known
within the energy resolution of the instrument.

The Principal Component Analysis method (PCA) has been used for the iterative track reconstruction algorithm. The first step is based on
the crystal signals only, without using the silicon detector information. Let $N$ be the total number of crystals with signal above a
given threshold, then the three-dimensional covariance matrix of the crystal coordinates $c_i^{(n)}$ (with $i=X,Y,Z$ and $n=1,\dots ,N$)
weighted with the signal amplitude is:
\begin{equation}
  M_{ij}^{CAL} = \frac{1}{\sum_{n=1}^N W^{(n)}_{CAL}} \sum_{n=1}^N W^{(n)}_{CAL} ( c_i^{(n)} - C_i ) ( c_j^{(n)} - C_j ) ,
\end{equation}
where $W^{(n)}_{CAL}$ is the assigned weight of the $n$-th crystal (equal to the signal $S_{CAL}^{(n)}$ on the $n$-th crystal) and $C_i$ is
the center of the distribution of the signals defined by:
\begin{equation}
  C_i = \frac{ \sum_{n=1}^N W^{(n)}_{CAL} c^{(n)}_i }{ \sum_{n=1}^N W^{(n)}_{CAL} }.
\end{equation}

As a first step, the eigenvector with the largest eigenvalue of the covariance matrix $M_{ij}^{CAL}$ is considered as estimator of the
primary track. The next step is to consider the energy deposits in the silicon detectors. At a given depth $Z$, a silicon micro-strip
detector measures only one coordinate ($X$ or $Y$), and to measure both typically two parallel detector planes arranged with perpendicular
segmentation directions are used. For this reason, the silicon measurements along the two segmentation directions are used to independently
reconstruct the two projections of the track on the $X-Z$ and $Y-Z$ views, and from these projections the track in the three-dimensional
space is then obtained. In the following the procedure for $X-Z$ is described, the one for $Y-Z$ being the same.

Let $J$ be the total number of strip with signal above threshold and with $d_i^{(k)}$  (with $i=X,Z$) denote the coordinate of the $k$-th
strip. The bidimensional covariance matrix for the PCA method is expressed by:
\begin{equation}
  M_{ij}^{SIL} = \frac{1}{\sum_{k=1}^J W^{(k)}_{SIL}} \sum_{k=1}^J W^{(k)}_{SIL} ( d_i^{(k)} - D_i ) ( d_j^{(k)} - D_j ),
\end{equation}
where $W^{(k)}_{SIL}$ is the weight for the $k$-th silicon readout channel, and with the centre of the distribution of the signals defined
by:
\begin{equation}
  D_i = \frac{ \sum_{k=1}^J W^{(k)}_{SIL} d^{(k)}_i }{ \sum_{k=1}^J W^{(k)}_{SIL} }.
\end{equation}

To improve the angular resolution, particular care must be taken in defining the weight $W^{(k)}_{SIL}$. In the algorithm used for this
work,
$W^{(k)}_{SIL}$ depends on the detector plane $P$ (with $P$ ranging from 1 to the number of silicon planes) to which the strip belongs.

Once a previous estimate of the particle track (in this case, the track obtained using only the crystal signals) is available, it can be
affirmed that the spatial dispersion for the $X$ coordinate of its intersection in the $P$-th plane can be approximated by a Gaussian
function. The standard deviation $\sigma_P$ can be estimated by simulation. The weighting factor is then defined by:
\begin{equation}
  W_{SIL}^{(k)} = \frac{ S^{(k)} }{ S_{P} } \times \frac{1}{\sqrt{2 \pi \sigma_P^2 }} \exp \left[ -\frac{1}{2} \left( \frac{d_i^{(k)} -
X_P}{\sigma_P} \right)^2 \right],
\end{equation}
where $X_P$ is the intersection of the track on the $P$-th plane given by the previous estimation, $S^{(k)}$ is the signal of the $k$-th
strip and $S_P$ is the total signal measured by the $P$-th silicon plane. The reason for the introduction of $S_P$ in the denominator is
that the angular dispersion in the development of the electromagnetic shower increases going deeper into the calorimeter due to multiple
scattering. Because the $S_P$ value increases as the amount of material crossed, the signals of silicon planes that lie deeper are
weighed to a lesser extent. The track reconstruction algorithm by means of the silicon detector signals uses the PCA method with the
covariance matrix $M_{ij}^{SIL}$. Since this new estimate depends on the previous one through the value of $W_{SIL}^{(k)}$, it can be
repeated by using the new estimate as the previous one until convergence.

\subsection{Results}
The angular resolution has been computed for simulated photons that convert into an e$^+$-e$^-$  pair in the first thin crystal, thus
producing a well-developed shower in the silicon detectors. The conversion efficiency in this case is $\sim$ 65\%. The resolution can be
improved by selecting only those events for which the tracking algorithm converged within a given number of iterative steps, reducing at
the same time the overall efficiency. This is an effective method to select ``golden samples'' of well-reconstructed events in offline
analysis, and for tuning the accuracy according to the needed total statistics. The resolution as a function of the overall efficiency
(conversion and selection) for 10 GeV photons is shown in figure \ref{fig:resvseff}. The TIC efficiency is the combined (conversion and
selection) efficiency, and tops at 65\% for 100\% selection efficiency. This can be regarded as an overestimation of the real
overall efficiency since other important factors are not considered in this estimate. For example, the trigger is expected to have a
sizeable impact, but for implementing a realistic trigger a more refined instrument design is needed, which would reasonably include other
sub-detectors (like an anticoincidence shield) which are unrelated to the problem of tracking gamma rays; so this aspect is neglected in
this study. For Fermi and DAMPE the efficiency has been estimated by dividing the effective areas quoted in \cite{FermiPerfWebPage} and
\cite{Chang2017} by the geometrical areas of the two instruments (a whole tracker plane for Fermi and the upper face of the BGO calorimeter
for DAMPE). At this energy the resolution of TIC is worse than Fermi (total) for the same efficiency of $\sim$ 36\% but can match Fermi by
reducing the efficiency to $\sim$ 14\%, and is generally better than DAMPE.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=14.5cm]{resvseff.eps}
\end{center}
\vspace{-0.5cm}
\caption{Angular resolution of TIC for 10 GeV photons as a function of the conversion+selection efficiency, compared to Fermi
\cite{FermiPerfWebPage} and DAMPE \cite{Chang2017}.}
\label{fig:resvseff}
\end{figure}

The performance of the tracking algorithm applied to the simulated data as a function of the gamma ray energy is compared to those of
Fermi and DAMPE in figure \ref{fig:rescomparison}. The general trend is a performance improvement with increasing energy, which is more
pronounced than that of Fermi and DAMPE and that provides considerably better performance above some ten GeV. In order to understand this
behavior it is useful to consider the limiting factors for the angular resolution at different energies. In the GeV region the tracking
precision is mainly limited by multiple scattering for both the pair conversion telescopes and TIC, the latter having thicker converters and
thus worse performance. As the energy increases the multiple scattering becomes less relevant and the pair conversion telescopes start to
get limited by the position resolution of their tracking systems, so that the tracking performance initially improves and then reaches a
plateau. On the contrary, the electromagnetic showers in TIC become more populated and less affected by statistical fluctuations, so that
the precision of the position measurements of the transverse profile continues to improve with increasing energy, and also does the tracking
precision. At 100 GeV the performance is better than Fermi and DAMPE, with comparable or better efficiency.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=14.5cm]{rescomparison.eps}
\end{center}
\vspace{-0.5cm}
\caption{Tracking performance of TIC according to Monte Carlo simulations, compared to Fermi \cite{FermiPerfWebPage} and DAMPE
\cite{Chang2017}.}
\label{fig:rescomparison}
\end{figure}

\section{Detector prototype}
\label{sect:prototype}

To validate the measurement principle and the simulation results, a small-scale prototype of the tracking calorimeter was built and
tested with particle beams. The idea behind this test was to compare the results obtained with real data with those obtained from a Monte
Carlo simulation of the prototype detector, in order to assess the validity of the full-scale simulation results described in section
\ref{sect:mcsimulations}. In order to reduce the construction costs, the detector was built by refurbishing the CaloCube prototype
\cite{Adriani2019} with additional layers of thin crystals obtained by sawing some spare CaloCube cubic scintillators, and with spare
one-sided silicon microstrip detectors from the DAMPE experiment \cite{Azzarello2016}.

The CaloCube crystals are made of CsI(Tl), which is a different material from the LYSO used for the full-scale simulations described in
section \ref{sect:mcsimulations}. For this reason, and also for the different sizes of the full-scale detector and the prototype, test beam
data have to be compared to a dedicated simulation of the prototype detector and not to the full-scale simulation. The dynamic range
required for the readout of the thin scintillators is narrower than that for cubic ones, since the former are less thick than the
latter and are placed on the external part of the instrument where the electromagnetic showers just start to develop; hence a lower amount
of scintillation photons with respect to the cubes is expected in thin elements. For this reason only the large-area Excelitas VTH2090
photodiode of the high-dynamic-range, dual-photodiode readout of CaloCube has been used for thin crystals. The DAMPE silicon sensors were
arranged in modules with two silicon wafers put side by side along the strip transversal direction; each wafer is 320 $\mu$m thick, with 384
readout strips (for a total of 768 strips per module) and a readout pitch of 240 $\mu$m. A small dead area of $\sim$2 mm was present between
the two active areas.

The prototype was designed to resemble the scheme in figure \ref{fig:design} with some adjustments due to experimental constraints, and is
schematically drawn in figure \ref{fig:prototype}. The {\it front} part is made of two trays of thin crystals with 5 elements along both X
and Y directions, for a total of 25 detectors per tray with 1.7 cm thickness, and three layers of silicon microstrip detectors with one
module on each layer. The {\it middle} part is made of three trays of 6x5 cubic elements, interleaved by two silicon layers. The front and
middle trays are placed ortogonally to the beam. The {\it back} part is made of 10 trays of cubes identical to those of the middle part but
arranged longitudinally with respect to the beam in two blocks of 5 trays each. All the silicon sensors have strips aligned along the Y
direction, so that only the shower development along X can be sampled with these detectors. Indeed, due to limited availability of spare
DAMPE sensors, it has not been possible to instrument also the Y direction. Thus, only the track projection on the XZ plane can be
reconstructed from the silicon data, and the spatially less-precise data from the scintillating crystals must be used for the YZ projection.
The whole structure was mounted on a rotating table which allowed to change the impact angle of the particle beam by rotating the detector
along the Y direction.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=14.5cm]{prototype.eps}
\end{center}
\vspace{-0.5cm}
\caption{Schematic view of the prototype as seen from top. The trays containing CsI(Tl) crystals are in blue, the silicon sensors in
orange. The beam impact angle could be changed by rotating the table on top of which the detector was mounted.}
\label{fig:prototype}
\end{figure}

\section{Test beam results}
\label{sect:tbresults}
The prototype was tested with electron beams of different energies at CERN test beam facilities (in May 2018 at PS and in September
2018 at SPS). The usage of electrons instead of gamma rays is not expected to introduce sizeable differences in the estimation of the
tracking resolution, since the tracking method is based on the sampling of electromagnetic showers which are very similar for electrons
and photons. Furthermore, it has sizeable experimental advantages, like the possibility to have monochromatic beams of known energy up to
$\sim$ 250 GeV and to use a very simple trigger system based on plastic scintillators read out by photomultiplier tubes.

Electron beams with energies ranging from 1 GeV to 100 GeV and incidence angles of 0 and 10 degrees were used during the two test sessions.
The rotating table supporting the detector was set up so that a nominal rotation angle of 0$^\circ$ corresponded to an alignment of the
detector axis with the nominal beam direction (i.e. to a nominal incidence angle of 0$^\circ$). Muons runs were used for calibrating the
response of the instrument and for aligning all the different detector elements.

The XZ projection of the track of the incident electrons was reconstructed with the algorithm described in section \ref{sect:algorithm}, and
compared with the true incidence angle estimated as the mean of all the reconstructed incidence angles at fixed positioning of the rotating
table. This procedure introduced a spread due to the different real impact angle of each electron with respect to the reference one, which
affected the estimated angular resolution, and that was taken into account with Monte Carlo simulations as described below.

An accurate simulation of the experimental setup was realized using FLUKA, including a detailed implementation of the detector geometry and
of several experimental effects. In the real setup the beam path was only partially in vacuum, with the electrons traveling through air
between disconnected segments of the beam pipe, and the traversed material (air, beam pipe windows and trigger scintillators) induced the
emission of forward-traveling bremsstrahlung photons which hit the detector and measurably affected the tracking performance. To account for
this effect, the beam line from the first in-air segment down to the detector was implemented in the simulation. Initial particle
kinematics was generated according to the angular and spatial distribution of the real beam as given by beam line simulations provided by
CERN. This allowed for reproducing the spread in the angular distributions described above. Relevant instrumental effects like the
capacitive couplings between the silicon microstrips (see \cite{Barbiellini2002} for a characterization of this effect for the silicon
sensors used in this study) were implemented as well.

The residual rotation misalignment of the silicon sensors might induce a spread in the measured angular resolution: since the Y coordinates
are always estimated using the scintillating crystals (due to the fact that the silicon sensors are segmented only along X) with a
resolution of about 1 cm, this big uncertainty in determining Y can affect the determination of X when paired with the residual tilt angle
of the silicon detectors. This effect was estimated by considering the size of the residual rotation given by the error on the
corresponding alignment parameter as obtained from a residuals minimization procedure. When paired with a Y resolution of 1 cm this
gives a spread of about 10 $\mu$m on X, which is small compared to the size of the residuals on X (about 50 $\mu$m at 100 GeV and 100
$\mu$m at 10 GeV), and thus the expected impact on the angular resolution of the track reconstructed with the silicon sensors is small. This
was verified with Monte Carlo simulations by applying a random shift on the X coordinates up to three times the expected value, with no
sizeable effect.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=7.5cm]{resolution1GeV.eps}
\includegraphics[width=7.5cm]{resolution5GeV.eps}
\end{center}
\vspace{-0.5cm}
\caption{Two-dimensional tracking resolution for 1 and 5 GeV electrons with normal incidence on the detector (red: Monte Carlo simulation,
black: test beam data).}
\label{fig:resolution1_5GeV}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=7.5cm]{resolution10GeV.eps}
\includegraphics[width=7.5cm]{resolution50GeV.eps}
\end{center}
\vspace{-0.5cm}
\caption{Two-dimensional tracking resolution for 10 and 50 GeV electrons with normal incidence on the detector (red: detector and upstream
simulation with instrumental effects and beam spread, blue: same as red with additional beam spread factor, black: test beam data).}
\label{fig:resolution10_50GeV}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=14.5cm]{resolution100GeV.eps}
\end{center}
\vspace{-0.5cm}
\caption{Two-dimensional tracking resolution for 100 GeV electrons with normal incidence on the detector (red: detector and upstream
simulation with instrumental effects and beam spread, blue: same as red with additional beam spread factor, black: test beam data).}
\label{fig:resolution100GeV}
\end{figure}

The agreement with test beam data is improved by accounting for these features, and actually it is reasonable that such a careful treatment
is needed to reproduce a very accurate angular resolution. The importance of such corrections becomes less relevant at lower energies where
the angular resolution becomes intrinsically worse and a very good agreement between test beam and Monte Carlo data can be obtained (figure
\ref{fig:resolution1_5GeV}). However, to fully reconcile test beam and Monte Carlo data at higher energies, an additional spread of about
0.04$^\circ$ had to be added to the nominal spread of the simulated beam; its contribution to the overall agreement becomes more and more
relevant as the energy increases (figures \ref{fig:resolution10_50GeV} and \ref{fig:resolution100GeV}). This correction accounts for
unknown residual systematic effects at all energies, but since its origin was not determined then it was not used when computing the
angular resolution from Monte Carlo data; the residual test beam - Monte Carlo resolution difference was taken as a systematic error.

The impact angle of the incident particle on the detector was found to have a negligible effect on the angular resolution. By comparing the
resolution plots for 0$^\circ$ and 10$^\circ$ incidence angles no sizeable difference was found; this was confirmed by Monte Carlo
simulations (see figure \ref{fig:resolution010deg}).

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=7.5cm]{resolution010degTB.eps}
\includegraphics[width=7.5cm]{resolution010degMC.eps}
\end{center}
\vspace{-0.5cm}
\caption{Two-dimensional tracking resolution for 100 GeV electrons for test beam (left) and Monte Carlo (right) data (black: 0$^\circ$,
blue: 10$^\circ$ incidence angle).}
\label{fig:resolution010deg}
\end{figure}

The angular resolution on the XZ projection plane are summarized in tables \ref{table:resolution68} and \ref{table:resolution95}. The
agreement between test beam and Monte Carlo data is quite satisfactory, especially at lower energies for 68\% containment of the
two-dimensional Point Spread Function (PSF). At 95\% PSF the agreement gets worse, likely due to residual instrumental effects that result
in a not-so-good modeling of the tails of the distribution by the Monte Carlo. Moreover, in the test beam analysis, no particular quality
selection criteria were applied in order reduce spurious events or the tails of the angular distribution.

\begin{table}
\begin{center}
\begin{tabular}[htbp]{|c|c|c|c|}
\hline
Energy (GeV) & Angle (deg)& Res. 68\% (deg) - TB & Res. 68\% (deg) - MC  \\
\hline
  1   &  0   & 3.72   $\pm$ 0.11    & 3.87 $\pm$ 0.12     \\
  5   &  0   & 0.985  $\pm$ 0.016   & 0.955 $\pm$ 0.014   \\
 10   &  0   & 0.4095 $\pm$ 0.0087  & 0.3971 $\pm$ 0.0050 \\
 50   &  0   & 0.1205 $\pm$ 0.0023  & 0.0995 $\pm$ 0.0050 \\
100   &  0   & 0.0897 $\pm$ 0.0010  & 0.0680 $\pm$ 0.0008 \\
100   & 10   & 0.0884 $\pm$ 0.0013  & 0.0646 $\pm$ 0.0004 \\
 \hline
\end{tabular}
\end{center}
\caption{Two-dimensional angular resolution in the XZ plane of the detector prototype, obtained from test beam (TB) and Monte
Carlo (MC) data, for 68\% PSF.}
\label{table:resolution68}
\end{table}
\begin{table}
\begin{center}
\begin{tabular}[htbp]{|c|c|c|c|c|c|}
\hline
Energy (GeV) & Angle (deg)& Res. 95\% (deg) - TB & Res. 95\% (deg) - MC  \\
\hline
  1   &  0   & 10.63 $\pm$ 0.15  & 11.44 $\pm$ 0.33 \\
  5   &  0   & 3.300 $\pm$ 0.063 & 3.282 $\pm$ 0.054 \\
 10   &  0   & 1.946 $\pm$ 0.062 & 1.790 $\pm$ 0.068 \\
 50   &  0   & 0.921 $\pm$ 0.029 & 0.710 $\pm$ 0.058 \\
100   &  0   & 0.678 $\pm$ 0.033 & 0.3077 $\pm$ 0.0099 \\
100   & 10   & 0.662 $\pm$ 0.036 & 0.259 $\pm$ 0.019 \\
 \hline
\end{tabular}
\end{center}
\caption{Two-dimensional angular resolution in the XZ plane of the detector prototype, obtained from test beam (TB) and Monte
Carlo (MC) data, for 95\% PSF.}
\label{table:resolution95}
\end{table}

\section{Conclusions}
In this paper, a new instrument design for tracking gamma rays in cosmic-ray space experiments has been proposed and characterized. The
TIC tracking calorimeter has comparable tracking performance with respect to the current generation of gamma-ray experiments (Fermi,
DAMPE) for energies from 1 GeV to $\sim$ 50 GeV, and better performance above $\sim$ 50 GeV, while allowing for important design
optimizations for multi-particle detection. Data acquired with a detector prototype at CERN test beam facilities showed a good agreement
with the Monte Carlo simulation, providing a deep insight on the relevant instrumental effects and a robust validation of the measurement
principle.

\acknowledgments
This work was financed by the INFN Commissione Scientifica Nazionale 5. The authors thank F. Cadoux for helping with realizing the detector,
and H. Wilkens, L. Gatignon, N. Charitonidis, B. Rae and M. Jeckel for test beam support.


\section{Introduction}
\IEEEPARstart{T}{opology} identification---which determines the connectivity of network nodes from an available measurement set---is vital in the electric grid, especially in distribution networks \cite{deka2017structure}, which may be switched between multiple operating configurations. Real-time topology awareness is critical for most control and optimization approaches and is essential for detecting unintentional changes caused by faulty equipment or cyberattacks. In three phase electrical networks, phase identification---determining the phase label of each phase at each bus---is part of topology identification. Correct phase labels are important for several applications, including allocating resources to minimize phase imbalance \cite{sun2015phase}. Together, topology and phase identification determine the entire network connectivity.
The literature demonstrates several approaches to phase identification---many heuristic---on simulated or real data sets. In \cite{arya2011phase} and \cite{jayadev2016novel} power balance constraints on load measurements are used to identify customer phase connectivity; this approach suffers from multiple feasible solutions. Another class of techniques applies clustering to power or voltage measurements to categorize phases, as in \cite{blakely2019spectral} and \cite{wang2016phase}. Correlations are a popular clustering distance metric for voltage magnitudes \cite{olivier2018phase, short2012advanced, pezeshki2012consumer}. In \cite{wen2015phase}, voltage correlations are used with voltage angle differences to match phases at two buses. Though some of these methods achieve good performance, they don't provide theoretical guarantees, stymieing understanding of why and how they work, and if and when they fail. For example, clustering algorithms such as k-means \cite{macqueen1967some} can suffer from local minima that may result in incorrect solutions \cite{jain1999data}. Furthermore, a theoretical justification enables understanding how load and network characteristics, such as the radial structure and line impedances, impact the success of an approach. \cite{wang2020maximum} presents a maximum likelihood estimator for phase identification, but it requires knowledge of the grid topology and is non-convex.

Topology identification has been extensively explored in the literature. Several works use voltage correlations or covariances to recover topology, based on heuristics \cite{bariya2018data}, theoretically justified patterns in covariances \cite{bolognani2013identification, deka2020graphical}, or conditional independence tests  \cite{weng2016distributed,deka2016estimating}. Another class of approaches uses linear regression to estimates the entire network model, including topology and line impedances, but require phasor current measurements in addition to voltages \cite{moffat2019unsupervised, yuan2016inverse, yu2017patopa, cavraro2018graph, park2020learning}. All these approaches assume that the system is balanced and can be approximated by a single phase network. This may lead to errors estimating unbalanced distribution networks. \cite{liao2019unbalanced, deka2019topology} present three-phase topology estimation, possibly with incorrect phase labels, using conditional independence or mutual information tests. However the sample requirements of these methods often exceed those of direct/greedy methods such as \cite{deka2017structure}. Post-topology recovery, \cite{liao2019unbalanced} uses a distance-based metric or mixed-integer program to classify unknown phases. This still falls short of a \emph{theoretically guaranteed}, \emph{polynomial-time algorithm} for phase recovery.

\textbf{Contribution:} This work presents a \textbf{g}reedy algorithm for joint \textbf{p}hase and \textbf{t}opology identification---termed \textbf{\textit{GPT}}---surpassing the prior work in that it
\begin{itemize}
    \item is provably correct under realistic assumptions and runs in polynomial time.
    \item is applicable to real networks, where some buses may have a subset of all three phases.
    \item requires \emph{only} voltage measurements; either phasors or magnitudes (with slight differences in performance).
    \item utilizes voltage \emph{statistics} enabling successful phase identification even from voltage magnitudes alone and in the presence of phase shifting transformers \cite{verboomen2005phase}.
\end{itemize}
GPT is based on a linearized, multi-phase power flow model mapping nodal current injections to nodal voltage phasors, and requires radial structure and diagonal dominance of line impedance matrices to guarantee correctness. When phases are known, GPT reduces to provable, greedy multi-phase topology learning generalizing prior work in \cite{deka2017structure} for the balanced setting. When topology is known apriori, it reduces to a local approach to phase identification. We demonstrate GPT's performance and improvement over prior work in both phase and topology recovery on multiple IEEE test networks simulated in Open-DSS, an open source distribution system simulator \cite{dugan2011open}.

The paper is organized as follows. Section \ref{sec:Model} presents a linearized model for unbalanced, three-phase networks with missing phases which is the theoretical basis of our work. Section \ref{sec:Phase} theoretically justifies a proximity metric for phase identification while Section \ref{sec:Topology} justifies a distance metric for topology identification once phases have been identified. Section \ref{sec:jointPhase_Top} puts the parts together to propose GPT: an algorithm for joint phase and topology recovery. Finally, Section \ref{sec:Results} presents validation on non-linear voltage data for three IEEE distribution test networks, simulated in OpenDSS. We compare performance with algorithms in prior work, demonstrating that our approach outperforms the prior work, and is robust to non-ideal measurements.
\begin{figure}
    \centering
    \includegraphics[width=.9\columnwidth]{figures/notation_toynet.png}
    \caption{(a) Notation visualized for two nodes and connecting line. (b) \emph{ToyNet}: A toy network with three, two, and one phase nodes used as a running example.}
    \label{fig:example}
\end{figure}

\subsection{Notation}
Notation is summarized below and visualized in Fig. \ref{fig:example}a.

\begin{table}[h!]
\begin{tabularx}{\columnwidth}{XX|XX}
 $\mathcal{N}$ & Node set  & $\mathcal{E}$ & Edge set  \\
 $\mathbf{V}_k = \begin{bmatrix} \mathbf{v}_k^a\\ \mathbf{v}_k^b\\ \mathbf{v}_k^c\end{bmatrix}$ & Node $k$ voltages & $\mathbf{I}_k = \begin{bmatrix} \mathbf{i}_k^a\\ \mathbf{i}_k^b\\ \mathbf{i}_k^c\end{bmatrix}$  &  Node $k$ current injections\\
 $\mathcal{M}_k$ & Node $k$ phases  & $\mathcal{M}_{kl}$  & Line $kl$ phases \\
 $\mathbf{V}$ & All nodal voltages & $\mathbf{I}$ & All nodal injections\\
 $\mathbf{Y}_{kl}$ & Multiphase admittance matrix  of line $kl$ & $\mathbb{Y}$ & Network admittance matrix\\
 $\mathbb{Y}_{k,l}$ & $\mathbb{Y}$ block for nodes $k$ \& $l$ &  $\mathbb{Y}_{k,l}^{\phi,\psi}$ & $\mathbb{Y}$ element for phase $\phi$ at node $k$ \& phase $\psi$ at node $l$\\
 $A^T$ & Transpose & $A^H$ & Conjugate transpose
\end{tabularx}
\end{table}

\section{Unbalanced Three Phase Model}\label{sec:Model}
Before introducing our unbalanced three phase network model, we review the \emph{single phase balanced} power flow model. In the single phase case, each line $ij$ has an associated scalar admittance $\mathbf{y}_{ij}$ and impedance $\mathbf{z}_{ij} = \mathbf{y}^{-1}_{ij}$. The single phase voltage phasors and current injections are related by the $(|\mathcal{N}| \times |\mathcal{N}|)$ system \emph{admittance matrix} $\mathbb{Y}$ as $\mathbf{I} = \mathbb{Y}\mathbf{V}$ with the form:
\begin{gather*}
    \mathbb{Y}_{i,j} = -\mathbf{y}_{ij},\ \mathbb{Y}_{i,i} = \sum_{ij \in \mathcal{E}} \mathbf{y}_{ij} \Rightarrow \mathbb{Y} = A\mathbf{D}A^T
\end{gather*}
$\mathbb{Y}$ can be factored into  $|\mathcal{N}|\times |\mathcal{E}|$ \emph{incidence} matrix $A$ and diagonal line impedance matrix $\mathbf{D}$ \cite{li1998laplacian}. Without loss of generality, we choose all edges to be directed away from the network ``root'', generally the point of common coupling (PCC) or substation. If edge $ij$ is oriented from $i$ to $j$, the corresponding elements of $A$ are
\begin{gather*}
    A_{i, ij} = 1,\ A_{j, ij} = -1, \ A_{r, ij} = 0 \text{~if~} r\neq i\neq j
\end{gather*}
where $A_{i, ij}$ is the $i^{th}$ element of the column corresponding to edge $ij$.
By definition, as $\mathbb{Y}1 = 0$, $\mathbb{Y}$ is not invertible. An invertible reduced admittance matrix, $\underline{\mathbb{Y}}$, is constructed by choosing a reference node $r$ and removing the corresponding row and column of $\mathbb{Y}$. Since the system is lossless, its inverse relates voltages and currents as follows:
\begin{align}
    \underline{\mathbb{Z}} = \underline{\mathbb{Y}}^{-1}, ~~ \underline{\mathbf{V}} = \underline{\mathbb{Z}}\underline{\mathbf{I}}
\end{align}
$\underline{\mathbf{V}}$  contains voltages differences to the reference voltage while $\underline{\mathbf{I}}$ contains current injections at non-reference nodes. Let $\mathcal{E}_i$ and $\mathcal{E}_j$ denote the edge sets on the unique path in the radial system to $r$ from nodes $i$ and $j$ respectively. The value of $\underline{\mathbb{Z}}_{i,j}$ is given by :
\begin{align*}
    \underline{\mathbb{Z}}_{i,j} = \sum_{kl \in (\mathcal{E}_i \cap \mathcal{E}_j)} \mathbf{z}_{kl}.
\end{align*}
Thus, the elements of $\underline{\mathbb{Z}}$ correspond to the impedances of common paths between node pairs and the reference \cite{moffat2019unsupervised,deka2017structure}.

The \emph{unbalanced three phase model} follows from the single phase one. To clarify definitions, we use \emph{ToyNet} (Fig. \ref{fig:example}b), a simple, unbalanced, three phase radial network, as a running example.
We begin with the model for a multiphase line $ij$, with phases $\mathcal{M}_{ij} \subseteq \{a,b,c\}$. The voltage across $ij$ is related to the current along each phase of the line by line impedance matrix $\mathbf{Y}_{ij}$:
\begin{align}\label{eq:line_impedance}
    \mathbf{I}_{ij} = \mathbf{Y}_{ij}(\mathbf{V}_i - \mathbf{V}_j)
\end{align}
$\mathbf{Y}_{ij}$ is the inverse of the {$(|\mathcal{M}_{ij}| \times |\mathcal{M}_{ij}|)$} line impedance matrix, $\mathbf{Z}_{ij}$. Eq. (\ref{eq:line_impedance}) for line $56$ in ToyNet is: {$\begin{bmatrix}\textbf{i}_{56}^a \\ \textbf{i}_{56}^c \end{bmatrix} = \mathbf{Y}_{56}\begin{bmatrix}\textbf{v}_5^a-\textbf{v}_6^a \\ \textbf{v}_5^c-\textbf{v}_6^c \end{bmatrix}$.} The node $i$ current injection, denoted $\mathbf{I}_i$, is a vector of injections on each phase of $i$, and is given by the sum of line flows: $\mathbf{I}_i = \sum_{ij \in \mathcal{E}} \mathbf{I}_{ij}$.
Building up from the current-voltage relations across individual lines in (\ref{eq:line_impedance}), the multi-phase voltages and currents injections across the network are related by $\mathbf{I} = \hat{\mathbb{Y}}\mathbf{V}$. Note this model can describe a network with a subset of phases at some nodes. $\hat{\mathbb{Y}}$ is the multi-phase system admittance matrix with dimensions $(\sum_{i \in \mathcal{N}} |\mathcal{M}_i|) \times (\sum_{i \in \mathcal{N}} |\mathcal{M}_i|)$. The $i,j$ \emph{block} of $\hat{\mathbb{Y}}$ is
\begin{align}\label{eq:3phaseYdef}
    \hat{\mathbb{Y}}_{i,j} = -\mathbf{Y}_{ij},\ \hat{\mathbb{Y}}_{i,i} = \sum_{ij \in \mathcal{E}} \mathbf{Y}_{ij}
\end{align}
Block $\hat{\mathbb{Y}}_{i,j}$ is $(|\mathcal{M}_i|\times |\mathcal{M}_j|)$, so $\mathbf{Y}_{ij}$ must be appropriately zero-padded or reduced if $i$ and $j$ don't have all the same phases. For ToyNet, $\hat{\mathbb{Y}}$ has the structure visualized in Fig. \ref{fig:modelViz}.
\begin{remark}
$\hat{\mathbb{Y}}$ can be factored into an incidence matrix $\hat{A}$, which captures the endpoints of each edge, and $\hat{\mathbf{D}}$, a block diagonal matrix of line admittances: $\hat{\mathbb{Y}} = \hat{A}\hat{\mathbf{D}}\hat{A}^T.$
$\hat{\mathbf{D}}$ has dimensions $(\sum_{ij \in \mathcal{E}}|\mathcal{M}_{ij}|) \times (\sum_{ij \in \mathcal{E}}|\mathcal{M}_{ij}|)$, with line admittance matrices $\mathbf{Y}_{ij}$ along the diagonal. $\hat{A}$ is $(\sum_{i \in \mathcal{N}}|\mathcal{M}_i|) \times (\sum_{ij \in \mathcal{E}}|\mathcal{M}_{ij}|)$ dimensional. Its rows correspond to phases at each bus, and columns to phases of each edge. With edges directed toward the root, assume edge $ij \in \mathcal{E}$ is oriented from $i$ to $j$. Then for every $ij \in \mathcal{E}$, with $\phi\in\mathcal{M}_{ij}$: $\hat{A}_{i, ij}^{\phi\phi} = 1,\ \hat{A}_{j, ij}^{\phi\phi} = -1$.
All other elements of $\hat{A}$ are zero. $\hat{A}$ for ToyNet is visualized in Fig. \ref{fig:modelViz}.
\end{remark}
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/ToyNet_Matrices_Small.png}
    \caption{Visualizing the structure of three phase admittance, adjacency and impedance matrices for ToyNet.}
    \label{fig:modelViz}
\end{figure}

\subsection{Inverting the model}
$\hat{\mathbb{Y}}$ maps voltages to current injections, but we use the inverse mapping for phase and topology identification. By definition (\ref{eq:3phaseYdef}), $\hat{\mathbb{Y}}$ is singular. Again, a reduction denoted $\underline{\hat{\mathbb{Y}}}$ is invertable. To obtain $\underline{\hat{\mathbb{Y}}}$, we remove the \emph{three} rows and columns of $\hat{\mathbb{Y}}$ corresponding to the three phases at reference node $r$. $\underline{\hat{\mathbb{Y}}}$ can be factored as: $\underline{\hat{\mathbb{Y}}} = \underline{A}\hat{\mathbf{D}}\underline{A}^T$ where $\underline{A}$ is obtained from $\hat{A}$ by deleting the three rows corresponding to $r$.
To derive $\underline{\hat{\mathbb{Y}}}^{-1}$, we begin with the right pseudoinverse of $\underline{A}^T$, which has the following properties.
\begin{lemma}\label{lemma:B}
Let $\underline{B}$ be the right pseudoinverse of $\underline{A}^T$, with rows corresponding to nodes and columns to edges. For $i \in \mathcal{N}$ with phases $\mathcal{M}_i$, let $\mathcal{E}_i$ be the edge set of the unique path to $r$. Then,
\begin{align*}
    \underline{B}_{i, kl}^{\phi\psi} =
\begin{cases} -1 ~~\forall \phi =\psi \in \mathcal{M}_i,\ \forall kl \in \mathcal{E}_i\\ 0 \textnormal{~~~~ otherwise } \end{cases}
\end{align*}
\end{lemma}
\begin{proof}
$\underline{B}_{i}^{\phi}$ is the row of $\underline{B}$ corresponding to phase $\phi$ at node $i$, while $\underline{B}_{,ij}^{,\phi}$ is the column corresponding to phase $\phi$ of edge $ij$. If $\phi \neq \psi$, $\underline{B}_{,ij}^{,\phi T}\underline{A}_{,kl}^{,\psi} =0$. A column of $\underline{A}$ has only two nonzero elements, so for $\phi=\psi$, we have
\begin{align*}
    \underline{B}_{,ij}^{,\phi T}\underline{A}_{,kl}^{,\phi} = -1\delta(ij\in\mathcal{E}_k) + 1\delta(ij\in\mathcal{E}_l)
\end{align*}
For any edge $kl \neq ij$, we will have either $(ij\in \mathcal{E}_k), (ij \in \mathcal{E}_l)$ \emph{or} $(ij\not\in \mathcal{E}_k), (ij \not\in \mathcal{E}_l)$. Thus $\underline{B}_{,ij}^{,\phi T}\underline{A}_{,kl}^{,\phi} =0$ for $kl \neq ij$. If $kl=ij$, we have $(ij\not\in \mathcal{E}_k), (ij \in \mathcal{E}_l)$ and $\underline{B}_{,ij}^{,\phi T}\underline{A}_{,ij}^{,\phi}=1$. Thus,
$\underline{B}_{,ij}^{,\phi T}\underline{A}_{,kl}^{,\psi} = 1 \textnormal{ iff }  ij=kl, \phi = \psi \Rightarrow \underline{B}^T\underline{A} = I.$

Now consider, $\underline{B}_{i}^{\phi}\underline{A}_{j}^{\psi T}$, the inner product of rows. If $\phi\neq \psi$, this is $0$. Consider when $\phi=\psi$ and $i\neq j$. If $j$ does not lie along the path from $i$ to the reference, then $jk \not\in \mathcal{E}_i$, i.e. there is no edge connected to $j$ in $\mathcal{E}_i$, and $\underline{B}_{i}^{\phi}\underline{A}_{j}^{\phi T} = 0$. In contrast, if $j$ lies along the path from $i$ to the reference, there \emph{must} be two edges $kj, jl \in \mathcal{E}_i$ oriented to and away from $j$ respectively, as the path passes through $j$. Then, $\underline{B}_{i}^{\phi}\underline{A}_{j}^{\phi T} = (-1\times -1) + (1\times -1) = 0$. If $i = j$, only edge $il =jl \in \mathcal{E}_i$, and  $\underline{B}_{i}^{\phi}\underline{A}_{i}^{\phi T} = 1$. Therefore, $\underline{B}_{i}^{\phi}\underline{A}_{j}^{\psi T} = 1 \textnormal{ iff } i=j, \phi=\psi$
\end{proof}
Therefore, ${\underline{\hat{\mathbb{Y}}}}^{-1} = \hat{\underline{\mathbb{Z}}}$ can be written as follows.
\begin{theorem}\label{thm:Y_inv}
The inverse of $\underline{\hat{\mathbb{Y}}}$ is given by:
\begin{align}\label{eq:Zdef}
    \underline{\hat{\mathbb{Z}}} = {\underline{\hat{\mathbb{Y}}}}^{-1} = \underline{B}\hat{\mathbf{D}}^{-1}\underline{B}^{T}
\end{align}
where $\hat{\mathbf{D}}^{-1}$ is block diagonal matrix of line impedance matrices: $\hat{\mathbf{D}}^{-1}_{ij,ij} = \mathbf{Z}_{ij}$. Further, the element of $\underline{\hat{\mathbb{Z}}}$ corresponding to phase $\phi$ at node $i$ and phase $\psi$ at node $j$ is given by:
\begin{align}\label{eq:Zelements}
    \underline{\hat{\mathbb{Z}}}_{ij}^{\phi\psi} = \sum_{kl \in (\mathcal{E}_i \cap \mathcal{E}_j)}\mathbf{Z}_{kl}^{\phi\psi}
\end{align}
\end{theorem}
\begin{proof}
Using the structure of $\underline{B}$ from Lemma \ref{lemma:B}, we have
$
    \underline{\hat{\mathbb{Z}}}\underline{\hat{\mathbb{Y}}} = (\underline{B}\hat{\mathbf{D}}^{-1}\underline{B}^{T})(\underline{A}\hat{\mathbf{D}}\underline{A}^{T}) = I
$. Thus $ \underline{\hat{\mathbb{Z}}} = {\underline{\hat{\mathbb{Y}}}}^{-1}$.

Consider the block of $\underline{\hat{\mathbb{Z}}}$ corresponding to nodes $i$ and $j$. Based on Lemma \ref{lemma:B}:
\begin{align}
    \underline{\hat{\mathbb{Z}}}_{ij}^{\phi\psi} &= \sum_{kl \in \mathcal{E}} \underline{B}_{i,kl}^{\phi\phi}\mathbf{Z}_{kl}^{\phi\psi}\underline{B}_{j,kl}^{\psi\psi T} &= \sum_{kl \in (\mathcal{E}_i \cap \mathcal{E}_j)}\mathbf{Z}_{kl}^{\phi\psi}\nonumber
\end{align}
\end{proof}
Intuitively, (\ref{eq:Zelements}) says that a change in current injection on phase $\psi$ at node $j$ will affect the voltage at phase $\phi$ at node $i$, proportional to the $(\phi\psi)$ impedance of the shared path $(\mathcal{E}_i \cap \mathcal{E}_j)$ from $i,j$ to $r$. In our definition, $\underline{\hat{\mathbb{Y}}}$ and $\underline{\hat{\mathbb{Z}}}$ are ordered with the phases of each node or edge grouped together. If all phases exist at all nodes, this is equivalent to a permutation of the three-phase model in \cite{deka2019topology}, where entries for one phase across all nodes and edges are grouped together.

Using this theorem, in the unbalanced three phase model, voltages are related to currents as:
\begin{align}\label{eq:3phaseOhms}
    \underline{\mathbf{V}} = \underline{\hat{\mathbb{Z}}}\underline{\mathbf{I}}
\end{align}
Here, $\underline{\mathbf{V}}$ contains nodal voltage \emph{differences} with the voltage for the \emph{matching phase} at the reference: $\underline{\mathbf{V}}_i^{\phi} = \mathbf{v}_i^{\phi}-\mathbf{v}_r^{\phi}$. $\underline{\mathbf{I}}$ contains all current injections except at the reference. Note that there is no assumption for all nodes having all three phases. In the next section, we will use the model of (\ref{eq:3phaseOhms}) to determine patterns in voltage statistics to enable phase and topology recovery.

\subsection{Current and Voltage Statistics}
We treat voltages as random variables driven by current via the model of (\ref{eq:3phaseOhms}). In our theoretical analysis, we assume
\begin{enumerate}[leftmargin=*]
\item current injections are uncorrelated across nodes and phases (including at a single node). \begin{align}\label{eq:assumption1}
        cov(\mathbf{i}_i^{\phi}, \mathbf{i}_j^{\psi}) \neq 0 \textnormal{ iff } (i=j) \cap (\phi=\psi)
    \end{align}
    \item current injections have equal variance at all nodes.
    \begin{align}\label{eq:assumption2}
        \forall i, \phi:\ var(\mathbf{i}_i^{\phi}) = s^2
    \end{align}
\end{enumerate}
As they are predominantly determined by loads---which are uncorrelated over time intervals on the order of seconds---current injections can be modeled as uncorrelated across nodes and phases when using high resolution measurements such as from PMUs. We assume PMUs report at $120 Hz$, but our methods apply if resolution is sufficient for measurements to be de-trended to remove inter-nodal correlations. Assumption (\ref{eq:assumption2}) is stronger but permissible in reasonably balanced networks. In Section \ref{sec:Results}, we evaluate how deviations from these assumptions impact recovery performance.

\section{Voltage Covariance for Phase Matching}\label{sec:Phase}
Voltages covariances are informative for phase identification. Under Assumptions (\ref{eq:assumption1},\ref{eq:assumption2}), the covariance of the voltage of phase $\phi$ at node $i$ and phase $\psi$ at node $j$, in the three-phase model (\ref{eq:3phaseOhms}), is given by:
\begin{align}\label{eq:covarv}
    cov(\mathbf{v}_i^{\phi}, \mathbf{v}_j^{\psi}) = cov(\underline{\hat{\mathbb{Z}}}_i^{\phi} \underline{\mathbf{I}},\underline{\hat{\mathbb{Z}}}_j^{\psi}\underline{\mathbf{I}}) = s^2Re((\underline{\hat{\mathbb{Z}}}_i^{\phi})^H\underline{\hat{\mathbb{Z}}}_j^{\psi})
\end{align}
We are interested in the sum of covariances for a particular \emph{phase ordering} between nodes $i$ and $j$. Consider the case where the phases at $i$ are a subset of those at $j$ ($\mathcal{M}_i \subseteq \mathcal{M}_j$). Let $\mathcal{O}$ denote the ordering/permutation of phases at $j$, where $\mathcal{O}(\phi)$ denotes the specific phase at $j$ matched to the phase $\phi$ at $i$. Then, the covariance sum for matching $\mathcal{O}$, denoted by $c_{ij}^{\mathcal{O}}$, is:
\begin{align}\label{eq:distMetric}
    c_{ij}^{\mathcal{O}} = \sum_{\phi \in \mathcal{P}_i} cov(\mathbf{v}_i^{\phi}, \mathbf{v}_j^{\mathcal{O}(\phi)})
\end{align}
Let $\underline{\hat{\mathbb{Z}}}_i$ denote the rows of $\underline{\hat{\mathbb{Z}}}$ corresponding to \emph{all} phases at node $i$, and $\underline{\hat{\mathbb{Z}}}_j^{\mathcal{O}}$ denote the rows corresponding to the phases at $j$ ordered according to $\mathcal{O}$. Then $c_{ij}^{\mathcal{O}}$ is:
\begin{align}
    c_{ij}^{\mathcal{O}} &= s^2Re(vec(\underline{\hat{\mathbb{Z}}}_i)^{H}vec(\underline{\hat{\mathbb{Z}}}_j^{\mathcal{O}})) = s^2Re(\sum_{k \in \mathcal{N}}vec(\underline{\hat{\mathbb{Z}}}_{ik}^H)vec(\underline{\hat{\mathbb{Z}}}_{jk}^{\mathcal{O}}))\nonumber\\
    &= s^2Re\Bigg[\sum_{k\in\mathcal{N}}\smashoperator[r]{\sum_{\substack{mn\in(\mathcal{E}_i \cap \mathcal{E}_k)\\pq\in(\mathcal{E}_j\cap\mathcal{E}_k)}}}vec(\mathbf{Z}_{mn}^{\mathcal{M}_i})^Hvec(\mathbf{Z}_{pq}^{\mathcal{O}})\Bigg]\label{eq:C_ijexpand}
\end{align}
The last equality follows from (\ref{eq:Zelements}). The contribution of a node $k$ to $c_{ij}^{\mathcal{O}}$ is the dot product of the common path lengths between $i,k$ and $j,k$. The following result shows how $c_{ij}^{\mathcal{O}}$ enables phase matching.

\begin{theorem}\label{thm:maxCov}
Consider $c_{ij}^{\mathcal{O}}$ given in (\ref{eq:distMetric}) for $\mathcal{M}_i \subseteq \mathcal{M}_j$. If condition (\ref{eq:lineCond}) holds for each pair of line impedance matrices, then $c_{ij}^{\mathcal{O}}$ is maximized when $\mathcal{O}$ corresponds to the correct phase matching between $i$ and $j$.
\begin{align}
    \forall \mathcal{M} \in \{\mathcal{M}_1,...,\mathcal{M}_n\},\ \forall st, kl \in \mathcal{E}:\nonumber\\ \mathcal{M} = \arg\max_{\mathcal{O}} Re\bigg[vec(\mathbf{Z}_{st}^{\mathcal{M}})^Hvec(\mathbf{Z}_{kl}^{\mathcal{O}(\mathcal{M})})\bigg]\label{eq:lineCond}
\end{align}
where $\mathcal{M}$ ranges over every nodal phase set $\mathcal{M}_i$.
\end{theorem}
\begin{proof} If (\ref{eq:lineCond}) holds, \emph{every} term in the summation in (\ref{eq:C_ijexpand}) is maximized by $\mathcal{O} = \mathcal{M}_i$. Therefore, $\mathcal{M}_i$ maximizes the sum, and $\mathcal{M}_i = \arg\max_{\mathcal{O}} c_{ij}^{\mathcal{O}}$
\end{proof}
(\ref{eq:lineCond}) is a condition on every pair of edges, $st$ and $kl$, in the network. It states that for every row subset of line impedance matrix $\mathbf{Z}_{st}$ (corresponding to each nodal phase set $\mathcal{M}_1,...,\mathcal{M}_n$), the matching rows of $\mathbf{Z}_{kl}$ produce the largest vectorized dot product. This is reasonable as real line impedance matrices are diagonally dominant. Condition (\ref{eq:lineCond}) depends on the particular network considered. In a network where all nodes have all three phases ($\mathcal{M}=\mathcal{M}_i = \{a,b,c\}$) the condition on the vectorized dot product involves the full three phase line impedance matrices. if some nodes have a subset of phases, it will involve sub-matrices of impedance matrices.
Note that, in general, $cov(\mathbf{v}_i^{\phi}, \mathbf{v}_j^{\psi}) > 0$ even if $\phi \neq \psi$. If node $i$ has phases $a,b$ and node $j$ has phases $a,c$ ($\mathcal{M}_i \not\subseteq \mathcal{M}_j$), minimizing $c_{ij}$ will incorrectly match $b, c$. GPT avoids such scenarios by ordering nodes; discussed later.
In summary, Theorem \ref{thm:maxCov} allows us to use $c_{ij}^{\mathcal{O}}$ as a proximity metric for phase matching.

\section{Voltage Difference Variances for Topology}\label{sec:Topology}
We use voltage difference variances for topology recovery. Define $d_{ij}$ to be the sum of the variance of the voltage differences between correctly matched phases of nodes $i,j$. Assuming $\mathcal{M}_i \subseteq \mathcal{M}_j$, $d_{ij}$ is:
\begin{align}\label{eq:sovDist}
     d_{ij} = \smashoperator[lr]{\sum_{\phi \in \mathcal{M}_i}}var(\mathbf{v}_i^{\phi}-\mathbf{v}_j^{\phi}) = \smashoperator[lr]{\sum_{\phi \in \mathcal{M}_i}}\mathbb{E}[(\mathbf{v}_i^{\phi}-\mathbf{v}_j^{\phi})-\mathbb{E}(\mathbf{v}_i^{\phi}-\mathbf{v}_j^{\phi})]^2
\end{align}
Lemma \ref{lem:minVarDiff} establishes trends in $d_{ij}$ along one phase.
\begin{lemma}\label{lem:minVarDiff}
Given the voltage on phase $\phi$ at node $i$:
\begin{align}\label{eq:minDist}
    \arg\min_{j} d_{ij}^{\phi} \triangleq
    \arg\min_{j} var(\mathbf{v}_i^{\phi} - \mathbf{v}_j^{\phi}) \in \textnormal{Parent/Child of } i
\end{align}
\end{lemma}
\begin{proof}
Expanding the difference, we obtain:
\begin{align}\label{eq:varsum}
    d_{ij}^{\phi} = var(\mathbf{v}_i^{\phi}-\mathbf{v}_j^{\phi}) = \sum_{n\in\mathcal{N}}\sum_{\psi\in\mathcal{M}_n}s_n^2 |\underline{\hat{\mathbb{Z}}}_{in}^{\phi\psi} - \underline{\hat{\mathbb{Z}}}_{jn}^{\phi\psi}|^2
\end{align}
where $s_n$ is the injection variance at node $n$. If the paths from nodes $k$ and $l$ to $r$ merge at node $n$, $\underline{\hat{\mathbb{Z}}}_{k,l}^{\phi,\psi} = \mathbf{e}_n^{\phi\psi}$, the impedance of the path from $n$ to $r$ along phase coupling $\phi, \psi$:
\begin{align}\label{eq:path}
\mathbf{e}_n^{\phi\psi} = \sum_{ij \in \mathcal{E}_n} \mathbf{Z}_{ij}^{\phi\psi}
\end{align}
To determine the minimizer of (\ref{eq:minDist}), consider two cases visualized in Fig. \ref{fig:AllRegions}. In case A, $j$ is the common ancestor of node $i,k$ on the path to the root. In case B, $j$ is an ancestor of $i$, while $k$ is an ancestor of $j$. In both cases, we show that $d_{ij}^{\phi},d_{jk}^{\phi} < d_{ik}^{\phi}$. Put together, for a given $i$, the minimizer $j$ of $d_{ij}^{\phi}$ is either the parent or child of $i$.

\textbf{Case A.} We split the sum in (\ref{eq:varsum}) into the regions $\mathcal{N}_i$ in Fig. \ref{fig:AllRegions}a. Using (\ref{eq:path}) in (\ref{eq:varsum}) for each region, we have
\begin{align*}
    d^{\phi}_{ik}-d_{ij}^{\phi} =
    \sum_{n \in \mathcal{N}_1,\psi \in \mathcal{M}_n} 0
    + \sum_{n \in \mathcal{N}_2,\psi \in \mathcal{M}_n}0 + \sum_{n \in \mathcal{N}_4,\psi \in \mathcal{M}_n}0\\
    + \sum_{n \in \mathcal{N}_3,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_j^{\phi\psi}-\mathbf{e}_n^{\phi\psi}|^2-|\mathbf{e}_j^{\phi\psi}-\mathbf{e}_j^{\phi\psi}|^2\bigg)\\
    + \sum_{n \in \mathcal{N}_5,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_j^{\phi\psi}-\mathbf{e}_k^{\phi\psi}|^2-|\mathbf{e}_j^{\phi\psi}-\mathbf{e}_j^{\phi\psi}|^2\bigg) > 0
\end{align*}
A similar argument shows $ d_{ik}^{\phi}-d_{kj}^{\phi}> 0$.\\
\textbf{Case B.} Now we split (\ref{eq:varsum}) over the regions in Fig. \ref{fig:AllRegions}b. Using (\ref{eq:path}), we have
\begin{align*}
    d_{ik}^{\phi}-d_{ij}^{\phi} =
    \sum_{n \in \mathcal{N}_1,\psi \in \mathcal{M}_n}0
    + \sum_{n \in \mathcal{N}_2,\psi \in \mathcal{M}_n}0\\
    + \sum_{n \in \mathcal{N}_3,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_n^{\phi\psi}-\mathbf{e}_k^{\phi\psi}|^2-|\mathbf{e}_n^{\phi\psi}-\mathbf{e}_n^{\phi\psi}|^2\bigg)\\
    + \sum_{n \in \mathcal{N}_4,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_j^{\phi\psi}-\mathbf{e}_k^{\phi\psi}|^2 - |\mathbf{e}_j^{\phi\psi}-\mathbf{e}_j^{\phi\psi}|^2\bigg)\\
    + \sum_{n \in \mathcal{N}_5,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_n^{\phi\psi}-\mathbf{e}_k^{\phi\psi}|^2-|\mathbf{e}_n^{\phi\psi}-\mathbf{e}_j^{\phi\psi}|^2\bigg)\\
    + \sum_{n \in \mathcal{N}_6,\psi \in \mathcal{M}_n}s_n^2\bigg(|\mathbf{e}_i^{\phi\psi}-\mathbf{e}_k^{\phi\psi}|^2-|\mathbf{e}_i^{\phi\psi}-\mathbf{e}_j^{\phi\psi}|^2\bigg) > 0
\end{align*}
A similar analysis shows $d_{ik}^{\phi}-d_{kj}^{\phi}> 0$. Thus the minimum is given by the parent/child of $i$.
\end{proof}
\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/allRegions.png}
    \caption{Regions of the radial network when (a) Case A: $k$ lies off the path between $i$ and the reference. (b) Case B: $i$ lies along the path from $j$ to the reference.}
    \label{fig:AllRegions}
\end{figure}Applying Lemma \ref{lem:minVarDiff} to all matched phases between two nodes gives the following result:
\begin{theorem}\label{thm:topology}
Given node $i$, the node $j$ which minimizes $d_{ij}$ in (\ref{eq:sovDist}) is either a parent or child of $i$.
\end{theorem}
If phases at each node are known, Thm. \ref{thm:topology} enables correct topology recovery with a greedy algorithm based on distance $d_{ij}$. Note that Lemma \ref{lem:minVarDiff} and Thm. \ref{thm:topology} hold for all uncorrelated injections even with unequal variances. Thus, Assumption \ref{eq:assumption2} can be relaxed for topology learning. We now have the tools for joint phase and topology recovery, detailed in the next section.

\section{Joint Phase \& Topology Identification}\label{sec:jointPhase_Top}
We propose GPT (Alg. \ref{alg:greedy}): a \emph{greedy} algorithm for joint phase and topology identification based on the nodal voltage properties of Sections \ref{sec:Phase}-\ref{sec:Topology}. GPT computes $c_{ij}^{\mathcal{O}}$'s (\ref{eq:covarv}) exhaustively (for all matching options), selecting maxima for phase matching (Theorem \ref{thm:maxCov}). Based on phase matchings, it computes $d_{ij}$'s (\ref{eq:sovDist}) exhaustively (for all node pairs), selecting minima for topology recovery (Theorem \ref{thm:topology}). GPT greedily builds a tree with node set $\mathcal{T}$, starting from node $i_0$ and iterating till all nodes have been added. In each iteration, a new node is added to the tree by choosing node $i \not\in \mathcal{T}$, which has the minimum value of $d_{ij}$ for all $j \in \mathcal{T}$, using the \emph{getNext} algorithm. When node $i$ is added to the tree, it's phases are determined based on their matching to the phases of node $j$.
\begin{algorithm}
\caption{\emph{GPT} Greedy topology and phase recovery}\label{alg:greedy}
\textbf{Input:} Multi-phase voltage time series for all nodes in $V$.\\
\textbf{Output:} $E$: network edge set, $P$: phase ordering of each node. \\
\begin{algorithmic}[1]
\State $N3,N2,N1 \gets$ set of three, two, one phase nodes in $V$.
\ForAll{$i \neq j \in N$}
\State $\forall \mathcal{O}$ compute $c_{ij}^{\mathcal{O}}$
\State $M[i, j] \leftarrow \arg\max c_{ij}^{\mathcal{O}}$ in (\ref{eq:distMetric}) \Comment{Phase matching}
\State $D[i, j] \leftarrow$ $d_{ij}$ in (\ref{eq:sovDist}) with matched $M[i, j]$
\EndFor
\State $T  \leftarrow \{i_0\}, N3 \leftarrow N3 \setminus i_0$,
 $P[i_0] \leftarrow [a, b, c]$ \Comment{Add first three phase node to tree $T$ and set phases}
\While{$N3 \not= \emptyset$}\Comment{Connect 3 phase nodes}
 \State $i$, $j$ $\leftarrow$ getNext($D$, $T$, $N3$)
 \State $T \leftarrow T \cup i, N3 \leftarrow N3 \setminus i$, $E \leftarrow E \cup e_{ij}$, {$P[i] \leftarrow M[i, j]$}
 \EndWhile
 \While{$N2 \not= \emptyset$}\Comment{Connect 2 phase nodes}
 \State $i$, $j$ $\leftarrow$ getNext($D$, $T$, $N2$)
 \State $T \leftarrow T \cup i$, $N2 \leftarrow N2 \setminus i$, $E \leftarrow E \cup e_{ij}$, {$P[i] \leftarrow M[i, j]$}
 \EndWhile
 \While{$N1 \not= \emptyset$}\Comment{Connect 1 phase nodes}
 \State $i$, $j$ $\leftarrow$ getNext($D$, $T$, $N1$)
 \State $T \leftarrow T \cup i$, $N1 \leftarrow N1 \setminus i$, $E \leftarrow E \cup e_{ij}$,
 {$P[i] \leftarrow M[i, j]$}
 \EndWhile
\end{algorithmic}
\end{algorithm}
\begin{algorithm}\label{alg:getNext}
\caption{\emph{getNext}}
\textbf{Input:} Pairwise distances in $D$, Added nodes in $T$, Nodes to add in $N$\\
\textbf{Output:} $i \in N$ to be added, $j \in T$ connected to $i$. \\
\begin{algorithmic}[1]
\State $d_{ij} \leftarrow \infty$, $i \leftarrow$ None, $j \leftarrow$ None\;
\ForAll{$b \in T, a \in N$}
\If{$D[a, b] < d_{ij}$}
\State $d_{ij} \leftarrow D[a, b], i \leftarrow a, j \leftarrow b$
\EndIf
\EndFor
\end{algorithmic}
\end{algorithm}
GPT adds $3$ phase, then $2$ phase, then $1$ phase nodes to the tree. The initial node must be three phase, making the reference an intuitive choice. By adding nodes in this order, GPT implicitly enforces the crucial fact that number of phases never increases moving from the substation to the network ends (a single phase node is never the parent of a three phase node) and avoids issues that can arise when applying a naive greedy algorithm to a network with a variable number of phases at each node. For example, suppose we are recovering the topology of \emph{ToyNet}. All nodes have been added to $\mathcal{T}$ except $6$, $7$, and $8$. To recover the correct topology, we should connect node $6$ to $5$ first. Then nodes $7$ and $8$ will get connected to $6$ naturally, as $d_{76} < d_{75}$ and $d_{86} < d_{85}$. However, consider $d_{65}$ and $d_{75}$,
\begin{align*}
    d_{65} = var(\mathbf{v}_6^a-\mathbf{v}_5^a) + var(\mathbf{v}_6^c-\mathbf{v}_5^c),\quad
    d_{75} = var(\mathbf{v}_7^c-\mathbf{v}_5^c)
\end{align*}
We have no guarantee that $d_{65} < d_{75}$ due to the presence of additional phase variance in $d_{65}$ illustrating how an algorithm that doesn't order nodes by decreasing number of phases may return an incorrect topology, unlike GPT.

\subsection{Alternative Estimation Scenarios}
In the general scenario, GPT recovers both phase and topology from voltage measurements. Our theoretical results also establish estimation methods for restricted settings:
\subsubsection{Phase Identification with Topology Information}
If topology is known, phases can be identified by greedily matching adjacent nodes using $c_{ij}^{\mathcal{O}}$ (\ref{eq:covarv}) across edges $ij \in \mathcal{E}$. \cite{olivier2018phase, pezeshki2012consumer, arghandeh2015topology} similarly use the Pearson correlation coefficient of voltages as the distance, which is related to the covariance but not theoretically justified. GPT is highly local, unlike \cite{blakely2019spectral, wang2016phase, olivier2018phase, short2012advanced, pezeshki2012consumer} which cluster all nodal voltages to recover phase. K-means is a popular clustering algorithm choice \cite{wang2016phase}. However, even if the correct phase matching is the globally optimal solution of the k-means cost, the optimization is non-convex and may not converge to the global minima. Our greedy approach, however, is guaranteed to result in the optimal solution.

\subsubsection{Topology Estimation with Phase Information}
If phase labels are known, $d_{ij}$ (\ref{eq:sovDist}) can be directly minimized to recover topology and GPT reduces to greedy spanning tree learning generalizing prior work for the single phase case \cite{deka2017structure}. Compared to \cite{liao2019unbalanced, deka2019topology} that use conditional independence tests and need matrix inversions, GPT has improved sample performance, as demonstrated in Section \ref{sec:Results}.

\subsubsection{Estimation using voltage magnitudes only}
While GPT is based on nodal voltage phasors, it can also use voltage magnitudes $v_i^{\phi} = |\mathbf{v}_i^{\phi}|$. This is theoretically justified by linearizing (\ref{eq:line_impedance}) for line $kl$:
\begin{align*}
    \mathbf{I}_{kl} &= \mathbf{Y}_{kl}\begin{bmatrix}e^{j\theta^a}(v_k^ae^{j\underline{\theta}_k^a}-v_l^ae^{j\underline{\theta}_l^a}) \\ e^{j\theta^b}(v_k^be^{j\underline{\theta}_k^b}-v_l^be^{j\underline{\theta}_l^b}) \\ e^{j\theta^c}(v_k^ce^{j\underline{\theta}_k^c}-v_l^ce^{j\underline{\theta}_l^c})  \end{bmatrix}\nonumber\\
    &\approx \mathbf{Y}_{kl}\mathbf{D}_r((V_k-V_l)+j(\underline{\theta}_k-\underline{\theta}_l))
\end{align*}
where $\theta^{\phi}$ is the phase $\phi$ reference angle and $\underline{\theta}_k^{\phi} = \theta_k^{\phi} - \theta^{\phi}$. The linearization assumes small magnitude deviations from the reference and small angle difference between neighboring nodes. Properties of voltage magnitudes across the network can then be derived under assumptions on  $\mathbf{Y}_{kl}\mathbf{D}_r, \mathbf{I}_{kl}$ to obtain GPT for phase and topology recovery.

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{figures/results/top_error_horiz.png}
\caption{GPT topology error vs number of samples for three test feeders and three noise levels. Samples are assumed to arrive at $120 Hz$. Insets show estimated (grey) and true (red) network lines across trials, with the opacity of grey lines indicating how many times the edge was recovered. GPT is evaluated on voltage phasor and magnitude data, with performance on magnitudes compared to state of the art in Liao2019 \cite{liao2019unbalanced}.}
\label{fig:top_results}
\end{figure*}
\begin{figure}[ht]
\centering
\includegraphics[width=\columnwidth]{figures/results/phase_error_horiz.png}
\caption{GPT phase matching error for three test feeders and three noise levels using voltage phasor and magnitude data. The last table allows comparison with the state of the art phase matching method of Olivier2018 \cite{olivier2018phase}. Mostly, GPT has comparable performance to Olivier2018, but significantly outperforms it on the 13 bus system.}
\label{fig:phase_results}
\end{figure}
\begin{figure}[ht]
\centering
\includegraphics[width=\columnwidth]{figures/results/corr_error_horz.png}
\caption{GPT topology error from voltage magnitudes as injections become increasingly correlated ($\epsilon \rightarrow 1$) for three test feeders and two noise levels. Notice the rapid rise in error with increasing correlation particularly for the largest network.}
\label{fig:corr_results}
\end{figure}
\section{Simulation Experiments}\label{sec:Results}
We present simulation results of GPT. We measure average errors in phase and topology recovery, normalized by network size:
\begin{align*}
    &\textnormal{Topology Error} = \frac{\textnormal{wrong edges} + \textnormal{missing edges}}{\textnormal{total edges}},\\
    &\textnormal{Phase Error} = \frac{\textnormal{wrong nodal phases}}{\textnormal{total nodal phases}}
\end{align*}
Further, we evaluate the algorithm's sensitivity to the following parameters.
\begin{itemize}[leftmargin=*]
    \item \emph{Measurement noise}: We add white noise $\mathbf{n}$ to original measurement $\mathbf{v}_i$: $\tilde{\mathbf{v}}_i= \mathbf{v}_i + \mathbf{n}$, defining $
        \textnormal{noise level}(\tilde{\mathbf{v}}_i) = \frac{var(\mathbf{n})}{var(\mathbf{v}_i)}
    $. As GPT uses voltage covariances, it depends on relative precision and not absolute accuracy and is immune to the stable transducer errors that afflict distribution PMU data \cite{brady2016uses}.
    \item \emph{Number of measurement samples}: Assuming $120$ Hz distribution PMU measurements \cite{von2014micro}, we record performance on $1$ second to $1$ minute of voltage data.
    \item \emph{Load Correlations}: We test GPT's sensitivity to the assumption of uncorrelated injections, by varying the correlations of the loads while maintaining their variance. This is done by setting load covariance matrix $\Sigma = \sigma^2 ((1-\epsilon) \mathbb{I} + \epsilon11^T)$. As $\epsilon \rightarrow 1$, injections become more correlated.
\end{itemize}

Three IEEE distribution test networks are simulated in OpenDSS: the 13 and 34 bus networks have some one and two phases buses, while the 37 bus network has all three phase buses \cite{kersting1991radial}. We modify the models by adding loads at every bus, and by disabling voltage regulators, which invalidate the assumption of voltages driven by injections. We fluctuate the load injections at each phase at each bus, and simulate the network with power flow to obtain non-linear voltages.

Fig. \ref{fig:top_results} plots topology recovery accuracy for three noise levels ranging from $0$ (no noise) to $10$, with $1$ second to $1$ minute of voltage magnitude or phasor measurements. PMUs are highly precise; and the noise level would realistically be $\sim 0.001$ \cite{brown2016characterizing, bariya2019empirical}. Nevertheless, GPT performs well under more noise as measurement samples increase. For all test networks and measurement durations, GPT achieves perfect topology recovery from voltage magnitudes for $0$ and $0.001$ noise. Insets in Fig. \ref{fig:top_results} show recovered topologies across trials. Note how errors are localized to a few nodes, and lower for voltage magnitudes. Fig \ref{fig:top_results} also compares performance on voltage magnitudes to Liao2019 \cite{liao2019unbalanced}, showing that GPT outperforms it across scenarios.

Fig. \ref{fig:phase_results} presents GPT's phase matching error on the same three networks and noise levels averaged across several sample durations (we found phase matching error to be invariant to sample duration). The error is compared to that of the approach in Olivier2018 \cite{olivier2018phase}. The methods have comparable performance, except on the 13 bus network, where GPT outperforms Olivier2018 across SNRs.

Fig. \ref{fig:corr_results} shows topology recovery sensitivity as injections stray from the uncorrelated assumption (\ref{eq:assumption1}). Error increases rapidly as loads become more correlated. In reality, over short time durations, it is reasonable to assume that injections will be uncorrelated across nodes or can be de-trended \cite{liao2019unbalanced}. We use at most one minute of data to recover phase and topology: short enough that the uncorrelated assumption should hold well.

Our polynomial time algorithms are suitable for real time application, taking on the order of seconds to recover phase and topology for the IEEE test cases shared here. On the largest 37 bus test case, the algorithm completes in $15$ seconds.

\section{Conclusion}
We presented GPT for joint phase and topology identification from voltage measurements in unbalanced three phase networks where each bus can have single, two, or three phases. GPT is polynomial time and relies on analytical trends in nodal voltage statistics. We proved its theoretical correctness under reasonable assumptions on load statistics and line impedances, demonstrating its efficacy on non-linear voltages from three test feeders simulated under realistic conditions. GPT is robust to measurement non-idealities, and outperforms the prior work in both phase and topology recovery.

Extensions of the work include handling of voltage regulators at some nodes and improving performance under load correlations. The linearized unbalanced three phase model presented here can motivate many future analytics, including extensions to settings with limited measurements, where our derived distance metrics could give a sense of the relative proximity of dispersed measurement points if not the precise topology. It is challenging to directly extend this work to non-radial settings as inter-nodal distances become difficult to define precisely in meshed networks. However, the model could inform graph theoretic monitoring approaches of meshed networks, similar to those in \cite{bariya2020physically}.

\bibliographystyle{ieeetr}