Great circle in UTM coordinates
The UTM coordinate system, is a kind of Transverse Mercator projection separated in longitude bands and restricted in latitud extent such that the distortions associated with the projection remain small. Also, the UTM coordinate system is conformal projection. Therefore, it preserves the angles. That means that within the UTM zones, a straight line in UTM projection is a good approximation great circle over the relatively short distances that could be found in a UTM zone. If you don't need great accuracy, you can just trace the straight line between P1 and P2.

The following plot shows the maximum deviation that a straight line have from a great circle. In this case P1 the red dot at the Equator and Longitude 3°E, and the maximum deviation is calculated for all points in the displayed box of 6°x6°.



While the differences remain small within distances of few hundred kilometers, bigger differences can be found over larger distances, specially if P1 and P2 are in different UTM zones.

If your application requires larger distances, you shouldn't use UTM projection. Alternatively, if you your application require greater accuracy over distances suitable for UTM projection, you can calculate the great circle trajectory, but I don't know any direct formulas.

What I have used, which I know is not what you are asking for, are some Matlab functions to convert to and from UTM coordinates and another that computes great circle tracks in geographic coordinates:

function [lat, lon] = greatCircle(P1,P2,N)
%GREATCIRCLE compute points along the great circle arc connecting points P1
%and P2
%   P1, P2 format: [lat,lon]
%   N: Number of points along the great circle arc

        lat1=P1(1);
        lon1=P1(2);

        lat2=P2(1);
        lon2=P2(2);

        lon12=lon2-lon1;

        a1=atan2d(cosd(lat2).*sind(lon12),(cosd(lat1)*sind(lat2)-sind(lat1)*cosd(lat1).*cosd(lon12)));
        %a2=atan2d(cosd(lat1).*sind(lon12),(-cosd(lat2)*sind(lat1)+sind(lat2)*cosd(lat1).*cosd(lon12)));

        a0=atand(sind(a1).*cosd(lat1)./sqrt((cosd(a1).^2)+(sind(a1).^2).*(sind(lat1).^2)));

        sigma01=atan2d(tand(lat1),cosd(a1));
        sigma12=atan2d(sqrt((cosd(lat1)*sind(lat2)-sind(lat1)*cosd(lat2)*cosd(lon12))^2+(cosd(lat2)*sind(lon12))^2),(sind(lat1)*sind(lat2)+cosd(lat1)*cosd(lat2)*cosd(lon12)));
        sigma02=sigma01+sigma12;

        lon01=atan2d(sind(a0)*sind(sigma01),cosd(sigma01));
        lon0=lon1-lon01;

        sigma=linspace(sigma01,sigma02,N);

        lat=atan2d(cosd(a0)*sind(sigma),sqrt(cosd(sigma).^2+sind(a0)^2*sind(sigma).^2));
        lon=wrapTo360(atan2d(sind(a0)*sind(sigma),cosd(sigma))+lon0);
end


And for the conversions:

function  [x,y,utmzone] = deg2utm(Lat,Lon,Huso)
% -------------------------------------------------------------------------
% [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
% Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.
%
% Inputs:
%    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
%    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
%
% Outputs:
%    x, y , utmzone.   See example
%
% Example 1:
%    Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
%    Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
%    [x,y,utmzone] = deg2utm(Lat,Lon);
%    fprintf('%7.0f ',x)
%       458731  407653  239027  230253  343898  362850
%    fprintf('%7.0f ',y)
%      4462881 5126290 4163083 3171843 4302285 2772478
%    utmzone =
%       30 T
%       32 T
%       11 S
%       28 R
%       15 S
%       51 R
%
% Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
%    LatDMS=[40 18 55.56; 46 17 2.04];
%    LonDMS=[-3 29  8.58;  7 48 4.44];
%    Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
%    Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
%    [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Author: 
%   Rafael Palacios
%   Universidad Pontificia Comillas
%   Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06, Aug/06
% Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern 
%    hemisphere coordinates. 
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------

% Argument checking
%
error(nargchk(2, 3, nargin));  %2 arguments required
CHuso=true;
if nargin==3
    CHuso=false;
end
n1=length(Lat);
n2=length(Lon);
if (n1~=n2)
   error('Lat and Lon vectors should have the same length');
end


% Memory pre-allocation
%
x=zeros(n1,1);
y=zeros(n1,1);
utmzone(n1,:)='60 X';

% Main Loop
%
for i=1:n1
   la=Lat(i);
   lo=Lon(i);

   sa = 6378137.000000 ; sb = 6356752.314245;

   %e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
   e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
   e2cuadrada = e2 ^ 2;
   c = ( sa ^ 2 ) / sb;
   %alpha = ( sa - sb ) / sa;             %f
   %ablandamiento = 1 / alpha;   % 1/f

   lat = la * ( pi / 180 );
   lon = lo * ( pi / 180 );
    if CHuso
        Huso = fix( ( lo / 6 ) + 31);
    end
   S = ( ( Huso * 6 ) - 183 );
   deltaS = lon - ( S * ( pi / 180 ) );

   if (la&lt;-72), Letra='C';
   elseif (la&lt;-64), Letra='D';
   elseif (la&lt;-56), Letra='E';
   elseif (la&lt;-48), Letra='F';
   elseif (la&lt;-40), Letra='G';
   elseif (la&lt;-32), Letra='H';
   elseif (la&lt;-24), Letra='J';
   elseif (la&lt;-16), Letra='K';
   elseif (la&lt;-8), Letra='L';
   elseif (la&lt;0), Letra='M';
   elseif (la&lt;8), Letra='N';
   elseif (la&lt;16), Letra='P';
   elseif (la&lt;24), Letra='Q';
   elseif (la&lt;32), Letra='R';
   elseif (la&lt;40), Letra='S';
   elseif (la&lt;48), Letra='T';
   elseif (la&lt;56), Letra='U';
   elseif (la&lt;64), Letra='V';
   elseif (la&lt;72), Letra='W';
   else Letra='X';
   end

   a = cos(lat) * sin(deltaS);
   epsilon = 0.5 * log( ( 1 +  a) / ( 1 - a ) );
   nu = atan( tan(lat) / cos(deltaS) ) - lat;
   v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
   ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2;
   a1 = sin( 2 * lat );
   a2 = a1 * ( cos(lat) ) ^ 2;
   j2 = lat + ( a1 / 2 );
   j4 = ( ( 3 * j2 ) + a2 ) / 4;
   j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
   alfa = ( 3 / 4 ) * e2cuadrada;
   beta = ( 5 / 3 ) * alfa ^ 2;
   gama = ( 35 / 27 ) * alfa ^ 3;
   Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
   xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;
   yy = nu * v * ( 1 + ta ) + Bm;

   if (yy&lt;0)
       yy=9999999+yy;
   end

   x(i)=xx;
   y(i)=yy;
   utmzone(i,:)=sprintf('%02d %c',Huso,Letra);
end


and

function  [Lat,Lon] = utm2deg(xx,yy,utmzone)
% -------------------------------------------------------------------------
% [Lat,Lon] = utm2deg(x,y,utmzone)
%
% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84).
% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez.
%
% Inputs:
%    x, y , utmzone.
%
% Outputs:
%    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
%    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
%
% Example 1:
% x=[ 458731;  407653;  239027;  230253;  343898;  362850];
% y=[4462881; 5126290; 4163083; 3171843; 4302285; 2772478];
% utmzone=['30 T'; '32 T'; '11 S'; '28 R'; '15 S'; '51 R'];
% [Lat, Lon]=utm2deg(x,y,utmzone);
% fprintf('%11.6f ',lat)
%    40.315430   46.283902   37.577834   28.645647   38.855552   25.061780
% fprintf('%11.6f ',lon)
%    -3.485713    7.801235 -119.955246  -17.759537  -94.799019  121.640266
%
% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds
% [Lat, Lon]=utm2deg(x,y,utmzone);
% LatDMS=dms2mat(deg2dms(Lat))
%LatDMS =
%    40.00         18.00         55.55
%    46.00         17.00          2.01
%    37.00         34.00         40.17
%    28.00         38.00         44.33
%    38.00         51.00         19.96
%    25.00          3.00         42.41
% LonDMS=dms2mat(deg2dms(Lon))
%LonDMS =
%    -3.00         29.00          8.61
%     7.00         48.00          4.40
%  -119.00         57.00         18.93
%   -17.00         45.00         34.33
%   -94.00         47.00         56.47
%   121.00         38.00         24.96
%
% Author: 
%   Rafael Palacios
%   Universidad Pontificia Comillas
%   Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------

% Argument checking
%
error(nargchk(3, 3, nargin)); %3 arguments required
n1=length(xx);
n2=length(yy);
n3=size(utmzone,1);
if (n1~=n2 || n1~=n3)
   error('x,y and utmzone vectors should have the same number or rows');
end
c=size(utmzone,2);
if (c~=4)
   error('utmzone should be a vector of strings like "30 T"');
end



% Memory pre-allocation
%
Lat=zeros(n1,1);
Lon=zeros(n1,1);


% Main Loop
%
for i=1:n1
   if (utmzone(i,4)&gt;'X' || utmzone(i,4)&lt;'C')
      fprintf('utm2deg: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n');
   end
   if (utmzone(i,4)&gt;'M')
      hemis='N';   % Northern hemisphere
   else
      hemis='S';
   end

   x=xx(i);
   y=yy(i);
   zone=str2double(utmzone(i,1:2));

   sa = 6378137.000000 ; sb = 6356752.314245;

%   e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
   e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
   e2cuadrada = e2 ^ 2;
   c = ( sa ^ 2 ) / sb;
%   alpha = ( sa - sb ) / sa;             %f
%   ablandamiento = 1 / alpha;   % 1/f

   X = x - 500000;

   if hemis == 'S' || hemis == 's'
       Y = y - 10000000;
   else
       Y = y;
   end

   S = ( ( zone * 6 ) - 183 ); 
   lat =  Y / ( 6366197.724 * 0.9996 );                                    
   v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
   a = X / v;
   a1 = sin( 2 * lat );
   a2 = a1 * ( cos(lat) ) ^ 2;
   j2 = lat + ( a1 / 2 );
   j4 = ( ( 3 * j2 ) + a2 ) / 4;
   j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
   alfa = ( 3 / 4 ) * e2cuadrada;
   beta = ( 5 / 3 ) * alfa ^ 2;
   gama = ( 35 / 27 ) * alfa ^ 3;
   Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
   b = ( Y - Bm ) / v;
   Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2;
   Eps = a * ( 1 - ( Epsi / 3 ) );
   nab = ( b * ( 1 - Epsi ) ) + lat;
   senoheps = ( exp(Eps) - exp(-Eps) ) / 2;
   Delt = atan(senoheps / (cos(nab) ) );
   TaO = atan(cos(Delt) * tan(nab));
   longitude = (Delt *(180 / pi ) ) + S;
   latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ...
                    (180 / pi);

   Lat(i)=latitude;
   Lon(i)=longitude;

end


I hope that helps.