List all degrees in decreasing order so that $94\ge d_1\ge d_2\ge\cdots\ge d_{95}$
For optimal graph, if there's edge (i,j) but no edge (j,k), (i,j,k are distinct vertexes) we have $d_i\gt d_k$, otherwise, replace edge (i,j) by edge (j,k) so could get another graph with larger degree square sum.
Now checking the adjacency matrix of the graph, we could find that for column i and column j of the matrix where $i\lt j$, except for elements in diagonal (i,i), if (j,x) is 1, (i,x) must be 1 too.
A valid candidate adjacency matrix looks like below
0 1 1 1 1 1 0 0
1 0 1 1 0 0 0 0
1 1 0 1 0 0 0 0
1 1 1 0 0 0 0 0
1 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0
Since we have $A\begin{pmatrix}1\\1\\\cdots\\1\end{pmatrix}=\begin{pmatrix}d_1\\d_2\\\cdots\\d_n\end{pmatrix}$ (where n=95), so we have
$\begin{pmatrix}1&1&\cdots&1\end{pmatrix} A^2\begin{pmatrix}1\\1\\\cdots\\1\end{pmatrix}=d_1^2+d_2^2+\cdots+d_n^2$ which is sum of all elements in $A^2$.
Now let's take a look at matrix $A^2$ and found element (i,j) of the matrix should be $min\{d_i,d_j\}$ or $min\{d_i,d_j\}-1$, and it is $min\{d_i,d_j\}-1$ iff the element (i,j) in A is 1, so we have
$d_1^2+d_2^2+...+d_n^2 = 1d_1+3d_2+5d_3+...+(2n-1)d_n -2e$, where e is total number of edges
or
$\begin{cases}n-1\ge d_1\ge d_2\ge\cdots\ge d_n\ge 0\\
d_1+d_2+\cdots+d_n=2e\\
d_1^2+d_2^2+\cdots+d_n^2 = d_1+3d_2+5d_3+\cdots+(2n-1)d_n-2e
\end{cases}$
More generally, some vertexes has same degrees we could merge them and generate new constrains
$\begin{cases}
0=i_0\lt i_1\lt i_2\lt\cdots\lt i_{t-1}\lt i_t=n\\
n-1\ge d_1\gt d_2\gt\cdots\gt d_{t-1}\gt d_t\ge d_{t+1}=0\\
(i_1-i_0)d_1+(i_2-i_1)d_2+\cdots+(i_t-i_{t-1})d_t=2e\\
(i_1-i_0)d_1^2+(i_2-i_1)d_2^2+\cdots+(i_t-i_{t-1})d_t^2= \\
(i_1^2-i_1-i_0^2+i_0)d_1+(i_2^2-i_2-i_1^2+i_1)d_2+\cdots+(i_t^2-i_t-i_{t-1}^2+i_{t-1})d_t
\end{cases}$
and find the maximal value of $(i_1-i_0)d_1^2+(i_2-i_1)d_2^2+\cdots+(i_t-i_{t-1})d_t^2$
Now let's allow $d_h$ to be real number since we only need get an upbound of it. Using Lagrange Multiplier, we could get
$
\begin{cases}
2d_h-u\times(i_{h-1}+i_h)+u-v=0&(2\le h\le t-1)\\
d_{h+1}+d_h-2u\times i_h+u-v=0&(1\le h\le t-1)
\end{cases}
$
and the range of h could include 1 or t if $d_1$ or $d_t$ does not take boundary value (n-1 or 0)
Let $e_h=d_h+\frac {u-v}2, j_h=u\times i_h$, we have
$\begin{cases}
2e_h=j_{h-1}+j_h&(2\le h\le t-1)\\
e_{h+1}+e_h=2j_h&(1\le h\le t-1)
\end{cases}$
and we could get $\begin{matrix}
2j_h=j_{h+1}+j_{h-1}&(2\le h\le t-2)
\end{matrix}$
More generally, if at least three $d_h$ does not take boundary value (n-1 or 0), for example $n-1\gt d_2 \gt d_3\gt d_4 \gt 0$ , and let's only treat $d_2,d_3,d_4,i_2,i_3$ as variable and all others as constant. We could use uppercase character to represent constant and lowercase variable to make the equations more readable, we have
$\begin{cases}
I_1\lt i_2 \lt i_3 \lt I_4 \\
D_1\gt d_2 \gt d_3 \gt d_4\gt D_5 \\
(i_2-I_1)d_2+(i_3-i_2)d_3 +(I_4-i_3)d_4 = A\\
(i_2-I_1)d_2^2+(i_3-i_2)d_3^2+(I_4-i_3)d_4^2 = (i_2^2-I_1^2)d_2+(i_3^2-i_2^2)d_3+(I_4^2-i_3^2)d_4+B
\end{cases}$
We need find maximal value of $(i_2-I_1)d_2^2+(i_3-i_2)d_3^2+(I_4-i_3)d_4^2= (i_2^2-I_1^2)d_2+(i_3^2-i_2^2)d_3+(I_4^2-i_3^2)d_4+B$
By Lagrange Multiplier, we have
$\begin{cases}
2(i_2-I_1)d_2-u(i_2^2-I_1^2)-v(i_2-I_1)=0\\
2(i_3-i_2)d_3-u(i_3^2-i_2^2)-v(i_3-i_2)=0\\
2(I_4-i_3)d_4-u(I_4^2-i_3^2)-v(I_4-i_3)=0\\
d_2^2-d_3^2-u(2i_2d_2-2i_2d_3)-v(d_2-d_3)=0\\
d_3^2-d_4^2-u(2i_3d_3-2i_3d_4)-v(d_3-d_4)=0
\end{cases}$
That's equivalent to
$\begin{cases}
2d_2-u(i_2+I_1)-v=0\\
2d_3-u(i_3+i_2)-v=0\\
2d_4-u(I_4+i_3)-v=0\\
d_2+d_3-2ui_2-v=0\\
d_3+d_4-2ui_3-v=0
\end{cases}$
And finally we get
$I_4-i_3=i_3-i_2=i_2-I_1=\frac{I_4-I_1}3=\Delta I\gt 0, d_2-d_3=d_3-d_4=\Delta d\gt 0$, where $\Delta I\gt 0$ is constant,and $\Delta d\gt 0$ is the a variable to be solved.
Let's use symbol $I_2=I_1+\Delta I, I_3=I_1+2\Delta I$, so where the local extreme value is reached, we have $i_2=I_2,i_3=I_3$
Next step we have $d_3=\frac{A}{3\Delta I}$ and
$(\frac{A}{3\Delta I}+\Delta d)^2+(\frac{A}{3\Delta I})^2+(\frac{A}{3\Delta I}-\Delta d)^2=(I_2+I_1)(\frac{A}{3\Delta I}+\Delta d)+(I_3+I_2)\frac{A}{3\Delta I}+(I_4+I_3)(\frac{A}{3\Delta I}-\Delta d)+\frac{B}{\Delta I}$
or $2\Delta d^2+4\Delta I \Delta d+H=0$, because $\Delta I\gt 0$,
If the equation of $\Delta d$ have real roots, at least one of them is negative. And if it has positive root, the negative one has larger absolute value.
And the result of target function becomes $\Delta I (d_2^2+d_3^2+d_4^2)=(\frac{A}{3\Delta I}+\Delta d)^2+(\frac{A}{3\Delta I})^2+(\frac{A}{3\Delta I}-\Delta d)^2=\frac{A^2}{3\Delta I^2}+2\Delta d^2$
So it takes larger value for negative root of $\Delta d$.
If we add constrain $i_2=I_2,i_3=I_3$ to original problem and ignore inequality constrain $D_1\gt d_2 \gt d_3 \gt d_4\gt D_5$ to get a problem
under constrain
$\begin{cases}
d_2+d_3 +d_4 = \frac{A}{\Delta I}\\
d_2^2+d_3^2+d_4^2 = (I_2+I_1)d_2+(I_3+I_2)d_3+(I_4+I_3)d_4+\frac{B}{\Delta I}
\end{cases}$
to get maximal value of $(I_2+I_1)d_2+(I_3+I_2)d_3+(I_4+I_3)d_4$
The geometry meaning of it is to find range of plane of given direction to intersect with a circle in 3D space. So there must be exact one maxima and one minima which is corresponding to the two roots of above equation of $\Delta d$ and the positive root must be corresponding to the minima. So it means in the original problem, the positive root of $\Delta d$ could not be the maxima too (but it could not e minima too). So it means the maximal value could only be reached in boundary condition (Smaller t required).
Next we need to enumerate situations with at most two $d_h$ doesn't take boundary value (n-1 and 0). It takes a lot of effort but we could verify all of them does not exceed 300000.
For example,
$n-1=D_1 \gt d_2 \gt d_3 \gt 0, 0=I_0\lt i_1 \lt i_2 \lt I_3=n$
We could result in $\begin{cases}D_1-d_2=d_2-d_3=\Delta d\\I_3-i_2=i_2-i_1=\Delta i\end{cases}$ and $\begin{cases}\Delta i\Delta d=\frac{n(n-1)-2e}3\\\Delta i^2-\frac35\Delta i-\frac{n(n-1)-2e}3=0\end{cases}$
The extreme value is $\frac53 (n(n-1)-2e)\Delta d-n(n-1)^2+4e(n-1)=246881.293\lt 300000$
And the global maxima is taken at
$n-1=D_1 \gt d_2 \gt 0, 0=I_0 \lt i_1 \lt I_2=n$,
we could get $d_2=\frac{2e-(n-1)i_1}{n-i_1}, i_1^2-(2n-1)i_1+2e=0$
and the global maxima value is $259781.511 \lt 300000$
Another approach is to use computer find exact solution:
Using the attribute of the adjacency matrix, we could easily get a formula to calculate the result efficiently by computer,
assume
f(n,E) is the max square sum of degrees of graph with n vertexes and E edges,
By removing the first column of and first row of the adjacency matrix, we could find that there're still $E-d_1$ edge left in a graph with only $a=d_1$ vertexs available so that $f(n,E)=\max_{\begin{cases}a(a+1)\ge 2E\\a\le n-1\end{cases}} \{f(a,E-a)+4E-3a+a^2\}$.
Using a simple C code below we could find f(95,2021)=258618.
#include <stdio.h>
#define MAXN 100
#define MAXE (MAXN*(MAXN-1)/2)
int ds[MAXN+1][MAXE+1];
int main()
{
int n,e,a;
ds[2][1]=2;
for(n=3;n<=MAXN;n++){
for(e=1;e<=n*(n-1)/2;e++){
a=n-1;if(a>e)a=e;
for(;a*(a+1)>=2*e;a--){
int v=ds[a][e-a]+4*e-3*a+a*a;
if(v>=ds[n][e])ds[n][e]=v;
}
printf("ds[%d,%d]=%d\n",n,e,ds[n][e]);
}
}
return 0;
}
And for E small enough such as E<n, it is easy to show that the best result is achieved when all edges are from a same vertex.
According to the formula $f(n,E)=\max_{\begin{cases}a(a+1)\ge 2E\\a\le n-1\end{cases}} \{f(a,E-a)+4E-3a+a^2\}$, it is likely that
$f(n,E)=f(n-1,E-n+1)+4E-3(n-1)+(n-1)^2$ for some n,E in reasonable range. Computer calcuation shows it work for n=95, it holds for $94\le E \le 2206$ (but not for any (n,E))
It means we could construct a graph with 24 vertexes with degree 94, one with degree 24+41=65, and 41 vertex with degree 25 and others 29 vertex with degree 24 to reach the sum square of degree 258618.
And it could be found that for most combination of (n,E) there's only a unique solution to reach the optimal result while some have 2~6 different solutions. But only one combination have 6 solutions:
(9,18)
ds[9,18]=192
011111111 (quasi-star)
101111111
110111000
111000000
111000000
111000000
110000000
110000000
110000000
011111111
101111111
110110000
111010000
111100000
110000000
110000000
110000000
110000000
011111111
101111000
110111000
111011000
111101000
111110000
100000000
100000000
100000000
011111110
101111110
110111110
111000000
111000000
111000000
111000000
111000000
000000000
011111100
101111100
110111100
111011100
111100000
111100000
111100000
000000000
000000000
011111100
101111100 (quasi-complete)
110111100
111011000
111101000
111110000
111000000
000000000
000000000