Nekomata, 10 bytes

ᵑʰ→:ṁ7ᵃc∙Å

Attempt This Online!

-2 bytes thanks to @Unrelated String.

This requires Nekomata v0.7.0.0.

Takes input as z [a,b,c].

ᵑʰ→:ṁ7ᵃc∙Å
ᵑ               Apply the following function z times:
 ʰ→               Increment any element of the list [a,b,c]
   :            Duplicate
    ṁ           Minimum
     7          Push 7
      ᵃ         Apply the following function to both the minimum and 7:
       c          Prepend to the duplicated list
        ∙       Dot product
         Å      Find the largest possible value

ʰ (\onAny) is a new built-in particle in Nekomata v0.7.0.0. It applies a function to any element of a list non-deterministically.

Nekomata v0.6.0.0, 12 bytes

3$ŧĦ+:ṁ7ᵃc∙Å

-2 bytes thanks to @Unrelated String.

Attempt This Online!

3$ŧĦ+
3$ŧ     Find a z-tuple of [0,1,2]
   Ħ    Histogram; count the occurrences of each of 0, 1, 2
    +   Add that to [a,b,c]

The other part is the same as the previous solution.

Python 3, 101 100 bytes

f=lambda a,b,c,z:max(f(a+1,b,c,z-1),f(a,b+1,c,z-1),f(a,b,c+1,z-1))if z else a*a+b*b+c*c+7*min(a,b,c)

Try it online!

Ruby, 78 bytes

f=->*a,z{z<1?a.sum{|x|x*x}+7*a.min: (0..2).map{|x|c=*a;c[x]+=1;f[*c,z-1]}.max}

Try it online!

05AB1E, 18 bytes

2ÝIãε2Ý¢+W7*snO+}à

Two separated inputs: \$z\$ and \$[a,b,c]\$.

Try it online or verify all test cases.

Explanation:

   ã             # Cartesian product
  I              # with the first input z, to create all z-sized lists using
2Ý               # list [0,1,2]
    ε            # Map over each inner list:
     2Ý¢         #  Count for each [0,1,2] how many times it occurs
        +        #  Add each to the second (implicit) input at the same position
         W       #  Push its minimum (without popping)
          7*     #  Multiply it by 7
            s    #  Swap so the triplet is at the top again
             n   #  Square each
              O  #  Sum them together
               + #  Add the min*7
    }à           # After the map: pop and leave the maximum
                 # (which is output implicitly as result)

K (ngn/k), 36 bytes

Find the max in all base cases. It takes two inputs a,b,c as a vector and z.

{|/(7*&/'k)+/'k*k:x+/:(y=+/')#+!3#3}

Try it online!

Charcoal, 34 bytes

⊞υＥ³ＮＦＮＦ⮌υＦ³⊞υＥκ⁺μ⁼λνＩ⌈Ｅυ⁺ΣＸι²×⁷⌊ι

Try it online! Link is to verbose version of code. Explanation:

⊞υＥ³Ｎ

Input the numbers of each of the scientific symbols as a list and push them to the list of lists.

ＦＮ

Loop over the number of special actions.

Ｆ⮌υＦ³⊞υＥκ⁺μ⁼λν

Generate all possibilities of an extra scientific symbol and add them to the list of lists of symbols.

Ｉ⌈Ｅυ⁺ΣＸι²×⁷⌊ι

Calculate the score of each of the lists and take the maximum.

In the case of one special action, it seems to be optimal to add it to your largest symbol, except in the case where the there is a unique smallest symbol and the difference between the largest and smallest symbol is less than 4, in which case you should add it to your smallest symbol. As a table:

x+1 x+1 x 1 -> x+1 x+1 x+1
x+2 x+1 x 1 -> x+2 x+1 x+1
x+2 x+2 x 1 -> x+2 x+2 x+1
x+3 x+1 x 1 -> x+3 x+1 x+1
x+3 x+2 x 1 -> x+3 x+2 x+1
x+3 x+3 x 1 -> x+3 x+3 x+1

Otherwise, increment the largest value. In the case of two special actions, repeating the above is optimal except in two cases; when the largest symbol is exactly one more than the other two then you should add one each to the latter, while when the largest symbol is exactly two more than the middle symbol and three more than the smallest symbol in which case you should add both to your largest symbol. As a table:

x+1 x   x 2 -> x+1 x+1 x+1
x+1 x+1 x 2 -> x+2 x+1 x+1
x+2 x+1 x 2 -> x+3 x+1 x+1
x+2 x+2 x 2 -> x+2 x+2 x+2
x+3 x+2 x 2 -> x+3 x+2 x+2
x+3 x+3 x 2 -> x+3 x+3 x+2

Otherwise, add two to the largest value. Note that the exceptions mentioned above are the first row of the above table and also x+3 x+1 x which should become x+5 x+1 x instead of x+4 x+1 x+1.

Arturo, 86 bytes

f:$[a,z]->?z=0->+∑map a=>[&^2]7*min a->max map 3'i[f map 3'j->+?i=j->1->0a\[j-1]z-1]

Try it!

Takes [a, b, c] as a block and z separately.

Explanation

f:$[a,z]->         ; a function f taking a block a and number z
?z=0->             ; is z zero? then...
+                  ; add
∑                  ; sum of
map a=>[&^2]       ; items of a squared
7*min a            ; and 7 times min of a
->                 ; otherwise...
max                ; max of
map 3'i[]          ; map over 1..3 and call it i
f map 3'j->...z-1  ; call f on <map over 1..3 and call it j> and z-1
+?i=j->1->0        ; add 1 if i=j; otherwise 0
a\[j-1]            ; to a[j-1]

JavaScript (Node.js), 85 bytes

f=(a,b,c,z)=>z?Math.max(f(b,c,a,z-=.25),f(a+1,b,c,~~z)):a*a+b*b+c*c+7*Math.min(a,b,c)

Try it online!

f=(a,b,c,z)=>z?               // Ending condition
    Math.max(
        f(b,c,a,z-=.25),      // Rotate
        f(a+1,b,c,~~z)        // Equals to max(
                              //   f(a+1,b,c,z-1),
                              //   f(b+1,c,a,z-1),
                              //   f(c+1,a,b,z-1),
       // Repeats but who cares => f(a+1,b,c,z-1),
                              //   f(b,c,a,z-1))
    )
:a*a+b*b+c*c+7*Math.min(a,b,c)

JavaScript (ES6), 91 bytes

f=(a,b,c,z,r=2)=>z?Math.max(f(a+1,b,c,z-1),r&&f(b,c,a,z,r-1)):a*a+b*b+c*c+7*Math.min(a,b,c)

Try it online!

Vyxal, 18 bytes

3Þ□↔Ṡ$v+ƛ²∑ng7*+;G

Try it Online!

This is a little ugly due to vyxal's vectorisation being silly, but it works.

   ↔               # Combinations of length z of
3Þ□                # Rows of 3x3 identity matrix
    Ṡ              # Sum each, giving a list of points to increment
     $v+           # Add the input to each
        ƛ       ;  # Over each
         ²∑        # Sum of squares
               +   # Plus
           ng      # minimum
             7*    # times 7
                 G # Take the largest value

Jelly, 19 16 bytes

3ṗ‘ṛ¦ƒ€µṂ×7+ḋ`)Ṁ

Try it online!

Basically the formula as given, all the way down to not actually exploiting the bound on \$z\$. Might still be a shorter way to do that--the best I have with Ṭ is 3Ṭ€¤ṗ+ƒ€µṂ×7+ḋ`)Ṁ for 17.

3ṗ                  [1 .. 3] to the Cartesian power of z. ([[]] if 0)
      €             For each set of indices to increment,
  ‘                 increment [a, b, c]
   ṛ¦               at
     ƒ              each index in sequence.
       µ      )     For each choice:
        Ṃ×7         multiply the min by 7
           +ḋ`      and add that to the sum of squares.
               Ṁ    Take the maximum of those results.

R, 60 bytes

f=\(x,z)`if`(z,max(apply(x+diag(3),2,f,z-1)),x%*%x+7*min(x))

Attempt This Online!

Takes input a,b,c as one vector and z separately.