What's the non-linearity of a BJT ? (about harmonic balance algorithm)

I'm recently reading about circuit simulation tutorials and quite confused by the concept of "non-linearity".

It's said that when performing harmonic balance, the circuit should be divided into 2 parts, "linear" and "non-linear", and a BJT, is absolutely "non-linear", reference I read is here.

My question is, when we're doing manual calculation for a BJT with <b,c,e> denoting the base, collector, emitter nodes, the only difference from other devices is the current control [LATEX]i_c=\beta i_b[/LATEX], and using Nodal Anlysis it's not quite different from a "linear" device.

Should I treat that [LATEX]i_c=a_0 +a_1 i_b +a_2 {i_b}^2 +...+ a_n {i_b}^n[/LATEX] for a rigorous simulation ? Although I don't know how to get the {[LATEX]a_n[/LATEX]} array, just a guess.

Any help will be appreciated !

A single transistor produces non-linear distortion of about 40% at a fairly high output and when it has no negative feedback. Negative feedback in a transistor circuit reduces distortion (but it also reduces voltage gain).

An OPA134 audio opamp is made with transistors and when negative feedback is applied then the non-linear distortion at any level is only 0.00008% or less.