A simple mass-spring-damper system is pretty simple. The differential equation can be simply written below:
$$ F= m \frac{d^2x}{dt^2} + b \frac{dx}{dt} + k x(t) $$
which is transformed via LaPlace to:
$$ \frac{F}{s}=x(s)(ms^2+bs+k) $$
Finding the damping value that makes the system critically damped is simple:
$$ s = \frac{-b ±\sqrt{b^2-4mk}}{2m} $$
therefore the system is critically damped when:
$b=\sqrt{4mk}$
Changing the spring to a progressive spring, the origional formula changes to:
$$ F= m \frac{d^2x}{dt^2} + b \frac{dx}{dt} + k x(t)^2. $$
How does the Laplace transform change from this and how does that effect the damping coefficient for critical damping?
Not sure about critically damped or expression for Laplace tranform without convolution in $s$ domain. But here is a way to obtain an analytic solution for $x(t)$.
The equation you want to solve for progressive spring case is not linear. So $e^{sx}$ is no longer an eigen function. I suggest you try to assume $x(t) = \sum_i a_i (t-r)^i$ i.e., $x(t)$ is analytic and solve for the equation $0=m\frac{d^2x}{dt^2} + b \frac{dx}{dt} + k * x(t)^2$.
The equation is $$\sum_{i\geq 2} m i(i-1) \times a_i \times (t-r)^{i-2}+\sum_{i\geq 1} bi \times a_i \times (t-r)^{i-1} + \sum_{i \geq 0} k \times \left(\sum_{k=0}^i a_{i-k} a_k \right) \times(t-r)^i = 0$$
Comparing coefficients of $(t-r)^0$, this will give $2ma_2 + b a_1 + k a_0^2 = 0$.
Comparing coefficients of $(t-r)^1$, this will give $6ma_3 + 2b a_2 + k (a_0a_1+a_1a_0) = 0$.
Comparing coefficients of $(t-r)^i$, this will give $$(i+2)(i+1)m a_{i+2} + (i+1)b a_{i+1} + k \left(\sum_{k=0}^i a_{i-k} a_k \right) = 0$$.
Start by assigning some value for $a_1,a_0$ and back solve to get a solution for $x(t)$. Although it seems there exists a solution if you solve like this, Not sure if its a useful solution for your purpose.
Note that $a_0 = x(r)$ and $a_1 = x'(r)$ which are exactly boundary conditions. You can assume $r=0$ for a solution near $t=0$. Note that the solution has a region of convergence and works only in an open interval around $r$. If you want solution for general analytic $F$, then replace the RHS with $b_i$ where $F = \sum_i b_i (t-r)^i$. Technically the last point is the usual trouble. You can only solve for $F$ which is analytic in an open interval around $r$.