how to define integral over d(sqrt(x))

I stuck how I define a function like below. Especially I want to know how to define the yellow box in the equation below.

Below is my script.

C0 = 0.05

def lin_equation(t, Ds):
    
    def integrand1(td):
        return Cs_func(t - td)
    

    def integrand2(td):
        return Cs_func(td)

    r = r_func(t)
    
    integral1, _ = quad(integrand1, 0, t**0.5)
    integral2, _ = quad(integrand2, 0, t)
    
    term1 = 2 * (Ds / np.pi)**0.5 * (C0 * np.sqrt(t) - integral1)
    term2 = (Ds / r) * (C0 * t - integral2)
    
    return term1 + term2

And if I run the code below,

Ds = 1e-10

t_lin = np.linspace(min(t_data), max(t_data), 100)
lin_eq = lin_equation(t_lin, Ds)

It gives me the error like below

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[455], line 4
      1 Ds = 1e-10
      3 t_lin = np.linspace(min(t_data), max(t_data), 100)
----> 4 lin_eq = lin_equation(t_lin, Ds)
      6 print(lin_eq)

Cell In[454], line 12, in lin_equation(t, Ds)
      8     return Cs_func(td)
     10 r = r_func(t)
---> 12 integral1, _ = quad(integrand1, 0, t**0.5)
     13 integral2, _ = quad(integrand2, 0, t)
     15 term1 = 2 * (Ds / np.pi)**0.5 * (C0 * np.sqrt(t) - integral1)

File C:\ProgramData\anaconda3\Lib\site-packages\scipy\integrate\_quadpack_py.py:438, in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst, complex_func)
    435     args = (args,)
    437 # check the limits of integration: \int_a^b, expect a < b
--> 438 flip, a, b = b < a, min(a, b), max(a, b)
    440 if complex_func:
    441     def imfunc(x, *args):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

So I have two questions.

First, It seems liek my code above is just integral dt, not integral d(t^0.5). How can I elaborate this?

Second, How can I fix that error?

Thank you

====================== Thanks to Matt's advice, I defined a new function to fix the value error issue. It does not give me errors, but I am unsure whether I am doing right. Could you check the code below and advise me?

def lin_equation(t, Ds):
    # Function to integrate Cs(t - td) over sqrt(td)
    def integrand1(tdd):
        return Cs_func(t - tdd**2)*0.5*tdd**-0.5
    
    # Function to integrate Cs(td) over td
    def integrand2(td):
        return Cs_func(td)

    r = r_func(t)
    
    integral1, _ = quad(integrand1, 0, t**0.5, limit=100)
    integral2, _ = quad(integrand2, 0, t, limit=100)
    
    term1 = 2 * (Ds / np.pi)**0.5 * (C0 * np.sqrt(t) - integral1)
    #term1 = 0
    term2 = (Ds / r) * (C0 * t - integral2)
    #term2 = 0
    
    return term1 + term2


def scalar_lin_eq(t_array,Ds):
    results_lin = np.empty(len(t_array))
    for i, t in enumerate(t_array):
        results_lin[i] = lin_equation(t, Ds)
        
    return results_lin

If the above function is correct, a new issue will arise. I am trying to fit this into my data. In short, How I define the Cs(t) function might be the problem? If you have any comments, I would appreciate it. Below is my script

popt_Ds, pcov_Ds = optimize.curve_fit(scalar_lin_eq, t_data, Gamma, p0=[1e-10])

t_lin = np.linspace(min(t_data), max(t_data), 100)
lin_eq = scalar_lin_eq(t_lin, popt_Ds[0])

#print(lin_eq)

plt.scatter(t_data, Gamma, label='Gamma')
plt.plot(t_lin, lin_eq, label='lin_eq', color='red')
plt.xlabel('time (t)')
plt.ylabel('Gamma')
plt.xscale("log")
plt.legend()
plt.show()

According to the literature, the equation in the first picture must be the S shape, but It seems like my script gives me monotonic increasing function.

The Cs(t) function is a cubic interpolated function from my data which is filtered by Gaussian_filter like below. The reason of the filter is that if I just run the simple interpolation like linear or cubic, the quad function reach to the limit, to It does not give me the correct value. So I tried to eliminate the noise.

Cs_smd = gaussian_filter(Cs_val, sigma=5)
Cs_cubic_interp = interpolate.interp1d(t_data, Cs_smd, kind='cubic', fill_value="extrapolate")

def Cs_func(t):
    return Cs_cubic_interp(t)

First, It seems liek my code above is just integral dt, not integral d(t^0.5).

Since you have both t and t_d in your equation, I think you mean it's like your code is just the integral with differential d(t_d), not d(sqrt(t_d)). Yes. To correct this, use the chain rule:

Substitute this into your integrand, and now the differential is dt_d, which quad is designed to handle, and make sure your limits of integration are defined in terms of t_d, not sqrt(t_d). Consider reviewing the Riemann-Stieltjes Integral and/or using the Mathematics Stack Exchange for more detailed math questions.

Second, How can I fix that error?

As mentioned in the comments, scipy.integrate.quad expects scalar limits of integration. You can avoid the error by passing the limits of integration one at a time. That is, you cannot simply pass the array t_lin because quad is not vectorized to handle that; you'd need to loop over the elements. Alternatively, if you don't mind using a private function that may be removed from SciPy at any time, you can replace your call to scipy.integrate.quad with a call to scipy.integrate._tanhsinh._tanhsinh. See "Is there a way to parallelize scipy.integrate.quad over a set of values?" for an example.